Skip to content

GitLab

S
starblade
Commit 1385b975 authored by Jakob Knollmueller's avatar Jakob Knollmueller
Browse files
Options

remove stuff

parent e091563a
Hide whitespace changes
Inline Side-by-side
Showing with 33 additions and 229 deletions
demos/KL_demo.py deleted 100644 → 0
View file @ e091563a
# 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, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2017-2018 Max-Planck-Society
# Author: Jakob Knollmueller
#
# Starblade is being developed at the Max-Planck-Institut fuer Astrophysik
import numpy as np
from astropy.io import fits
from matplotlib import pyplot as plt
from multiprocessing import Pool
import nifty4 as ift
from nifty4.library.nonlinearities import PositiveTanh
import starblade as sb
from starblade.starblade_energy import StarbladeEnergy
from starblade.starblade_kl import StarbladeKL
def power_update(KL_energy):
power = 0.
for energy in KL_energy.energy_list:
power += ift.power_analyze(FFT.inverse_times(energy.s),
binbounds=p_space.binbounds)
power /= len(KL_energy.energy_list)
return power
if __name__ == '__main__':
#specifying location of the input file:
path = 'data/hst_05195_01_wfpc2_f702w_pc_sci.fits'
path = 'data/frame-u-006174-2-0094.fits'
# path = 'data/frame-g-002821-6-0141.fits'
path = 'data/frame-g-007812-6-0100.fits'
path = 'data/frame-i-004874-3-0692.fits'
# data = fits.open(path)[1].data
data = fits.open(path)[0].data#[1000:,1250:]
data -= data.min() - 0.001
# data = np.exp(2*(1.-plt.imread('data/sdss.png').T[0]))
# data = (plt.imread('data/m51_3.jpg').T[0])
# data = (plt.imread('data/12_FBP.png').T[0])
#
# data = data.clip(min=0.001)
data = np.ndarray.astype(data, float)
vmin = np.log(data.min()+0.01)
vmax = np.log(data.max())
plt.imsave('data.png', np.log(data))
postanh=PositiveTanh()
alpha = 1.5
s_space = ift.RGSpace(data.shape, distances=len(data.shape) * [1])
h_space = s_space.get_default_codomain()
data = ift.Field(s_space,val=data)
FFT = ift.FFTOperator(h_space, target=s_space)
binbounds = ift.PowerSpace.useful_binbounds(h_space, logarithmic = False)
p_space = ift.PowerSpace(h_space, binbounds=binbounds)
initial_spectrum = ift.power_analyze(FFT.inverse_times(ift.log(data)),
binbounds=p_space.binbounds)
initial_spectrum /= (p_space.k_lengths+1.)**4
update_power = True
initial_x = ift.Field(s_space, val=-1.)
alpha = ift.Field(s_space, val=alpha)
q = ift.Field(s_space, val=1e-30)
ICI = ift.GradientNormController(iteration_limit=100,
tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=ICI)
parameters = dict(data=data, power_spectrum=initial_spectrum,
alpha=alpha, q=q,
inverter=inverter, FFT=FFT,
newton_iterations=5, update_power=update_power)
current_x = initial_x
for i in range(10):
Starblade = StarbladeEnergy(position=current_x, parameters=parameters)
samples = []
for i in range(3):
sample = Starblade.curvature.inverse.draw_sample()
samples.append(sample)
problem = StarbladeKL(current_x, samples,parameters)
controller = ift.GradientNormController(name="Newton",
tol_abs_gradnorm=1e-5,
iteration_limit=5)
minimizer = ift.RelaxedNewton(controller=controller)
problem, convergence = minimizer(problem)
current_x = problem.position
parameters['power_spectrum'] = power_update(problem)
Starblade = StarbladeEnergy(position=current_x, parameters=parameters)
# Starblade = sb.build_starblade(data, alpha=alpha)
# for i in range(10):
# Starblade = sb.starblade_iteration(Starblade)
#
# #plotting on logarithmic scale
plt.imsave('diffuse_component.png', (Starblade.s).val,vmin=vmin, vmax=vmax)
plt.imsave('pointlike_component.png', Starblade.u.val, vmin=vmin, vmax=vmax)
Starblade = StarbladeEnergy(position=current_x, parameters=parameters)
var = 0.
mean = 0
samps = 30
for i in range(samps):
sam = postanh(Starblade.position+Starblade.curvature.inverse.draw_sample())
mean += sam
var += sam**2
var /= samps
mean /= samps
var -= mean**2
mask = ift.sqrt(var) < 0.01 +0.
plt.imsave('masked_points.png', mask.val * Starblade.u.val, vmin=vmin, vmax=vmax)
plt.imsave('masked_diffuse.png', mask.val * Starblade.s.val)
plt.imsave('std.png', np.log(np.sqrt(var.val)*data.val), vmin=-3.3)
# plt.figure()
# k_lenghts = Starblade.power_spectrum.domain[0].k_lengths
# plt.plot(k_lenghts, Starblade.power_spectrum.val)
# plt.title('power spectrum')
# plt.yscale('log')
# plt.xscale('log')
# plt.ylabel('power')
# plt.xscale('harmonic mode')
# plt.savefig('power_spectrum.png')
demos/clipping.py deleted 100644 → 0
View file @ e091563a
# 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, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2017-2018 Max-Planck-Society
# Author: Jakob Knollmueller
#
# Starblade is being developed at the Max-Planck-Institut fuer Astrophysik
import numpy as np
from astropy.io import fits
from matplotlib import pyplot as plt
from scipy.ndimage.filters import median_filter
import starblade as sb
if __name__ == '__main__':
#specifying location of the input file:
# path = 'data/hst_05195_01_wfpc2_f702w_pc_sci.fits'
# data = fits.open(path)[1].data
path = 'data/frame-i-004874-3-0692.fits'
path ='data/check.fits'
# data = fits.open(path)[1].data
data = fits.open(path)[0].data[1000:,1250:]
data -= data.min() - 0.001
data = data.clip(min=0.001)
data_true = data.copy()
data = np.ndarray.astype(data, float)
vmin = np.log(data.min()+0.01)
vmax = np.log(data.max())
local_size = 4
for i in range(5):
for i in range(data.shape[0]/local_size):
for j in range(data.shape[1]/local_size):
local_data = data[i*local_size:(1+i)*local_size,j*local_size:(1+j)*local_size]
local_data_median = np.median(local_data)
local_data_var = local_data.var()
local_data = local_data.clip(min=local_data_median - 3*np.sqrt(local_data_var),
max=local_data_median + 3*np.sqrt(local_data_var))
data[i * local_size:(1 + i) * local_size, j * local_size:(1 + j) * local_size] = local_data
background = np.empty_like(data)
crowded = np.zeros_like(data)
for i in range(data.shape[0] / local_size):
for j in range(data.shape[1] / local_size):
local_true_data = data_true[i * local_size:(1 + i) * local_size, j * local_size:(1 + j) * local_size]
local_data = data[i * local_size:(1 + i) * local_size, j * local_size:(1 + j) * local_size]
local_true_var = local_true_data.var()
local_var = local_data.var()
if 0.8 * np.sqrt(local_true_var) > np.sqrt(local_var):
background[i * local_size:(1 + i) * local_size,
j * local_size:(1 + j) * local_size] = 2.5*np.median(local_data)-1.5*local_data.mean()
crowded[i * local_size:(1 + i) * local_size,
j * local_size:(1 + j) * local_size] = 1.
else:
background[i * local_size:(1 + i) * local_size,
j * local_size:(1 + j) * local_size] = local_data.mean()
background = median_filter(background, size=(local_size,local_size))
# alpha = 1.25
# Starblade = sb.build_starblade(data, alpha=alpha)
# for i in range(10):
# Starblade = sb.starblade_iteration(Starblade)
#
# plotting on logarithmic scale
# background += background.min()
plt.gray()
plt.imsave('diffuse_component.png', np.log(background))#, vmin=vmin, vmax=vmax)
plt.imsave('pointlike_component.png', (data_true - background), vmin=vmin, vmax=vmax)
plt.imsave('crowded.png',crowded)
# plt.figure()
# k_lenghts = Starblade.power_spectrum.domain[0].k_lengths
# plt.plot(k_lenghts, Starblade.power_spectrum.val)
# plt.title('power spectrum')
# plt.yscale('log')
# plt.xscale('log')
# plt.ylabel('power')
# plt.xscale('harmonic mode')
# plt.savefig('power_spectrum.png')
starblade/starblade_kl.py
View file @ 1385b975
......@@ -20,6 +20,30 @@ from nifty4 import Energy, Field, DiagonalOperator, InversionEnabler
from starblade_energy import StarbladeEnergy
class StarbladeKL(Energy):
"""The Kullback-Leibler divergence for the starblade problem.
Parameters
----------
position : Field
The current position of the separation.
samples : List
A list containing residual samples.
parameters : Dictionary
Dictionary containing all relevant quantities for the inference,
data : Field
The image data.
alpha : Field
Slope parameter of the point-source prior
q : Field
Cutoff parameter of the point-source prior
power_spectrum : callable or Field
An object that contains the power spectrum of the diffuse component
as a function of the harmonic mode.
FFT : FFTOperator
An operator performing the Fourier transform
inverter : ConjugateGradient
the minimization strategy to use for operator inversion
"""
def __init__(self, position, samples, parameters):
super(StarbladeKL, self).__init__(position=position)
......
starblade/sugar.py
View file @ 1385b975
......@@ -80,6 +80,8 @@ def starblade_iteration(starblade, samples=3):
----------
starblade : StarbladeEnergy
An instance of an Starblade Energy
samples : int
Number of samples drawn in order to estimate the KL. If zero the MAP is calculated (default: 3).
"""
controller = ift.GradientNormController(name="Newton",
tol_abs_gradnorm=1e-8,
......@@ -152,6 +154,13 @@ def multi_starblade_iteration(MultiStarblade, processes = 1):
return NewStarblades
def update_power(energy):
""" Calculates a new estimate of the power spectrum given a StarbladeEnergy or StarbladeKL.
For Energy the MAP estimate of the power spectrum is calculated and for KL the variational estimate.
----------
energy : StarbladeEnergy or StarbladeKL
An instance of an StarbladeEnergy or StarbladeKL
"""
if isinstance(energy, StarbladeKL):
power = 0.
for en in energy.energy_list:
......
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