From 75f98c034a511acf1124543e27eaa7f195befdc9 Mon Sep 17 00:00:00 2001
From: "Joern Ungermann (IEK-7 FZ-Juelich)" <j.ungermann@fz-juelich.de>
Date: Thu, 30 Oct 2014 17:22:10 +0100
Subject: [PATCH] fixed norming of gradient for lp minimization.
---
pyramid/reconstruction.py | 2 +-
scripts/collaborations/example_joern.py | 59 ++++++++++++++++++++-----
2 files changed, 49 insertions(+), 12 deletions(-)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index b5a1500..4f85114 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -153,7 +153,7 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
lp = regularisator.norm
lq = jnorms.LPPow(q, 1e-20)
def preconditioner(_, direc):
- direc_p = direc / direc.max()
+ direc_p = direc / abs(direc).max()
direc_p = 10 * (1. / q) * lq.jac(direc_p)
return direc_p
diff --git a/scripts/collaborations/example_joern.py b/scripts/collaborations/example_joern.py
index 4723f25..20e089b 100644
--- a/scripts/collaborations/example_joern.py
+++ b/scripts/collaborations/example_joern.py
@@ -7,7 +7,8 @@ import os
import numpy as np
import pickle
-
+import matplotlib
+matplotlib.use("Agg")
import matplotlib.pyplot as plt
import pyramid
@@ -17,18 +18,44 @@ from pyramid.phasemapper import pm
from pyramid.dataset import DataSet
from pyramid.regularisator import ZeroOrderRegularisator, FirstOrderRegularisator
import pyramid.reconstruction as rc
-
from time import clock
import logging
import logging.config
+if "PBS_ARRAYID" in os.environ:
+ job_number = int(os.getenv("PBS_ARRAYID"))
+ experiments = []
+ for lam in [1e-9, 1e-8, 1e-7, 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1]:
+ for order in [0, 1]:
+ for p in [2, 1.5, 1.1, 1.01, 1.001]:
+ for use_mask in [True, False]:
+ experiments.append(
+ (lam, order, p, use_mask))
+ assert 0 <= job_number < len(experiments), \
+ "job id too large, maximum={}".format(len(experiments) - 1)
+ lam, order, p, use_mask = experiments[job_number]
+ print experiments[job_number]
+ dirname = "new-order_{}-p_{}-mask_{}-lambda_{:.9f}-job{}/".format(
+ order, p, use_mask, lam, job_number)
+ if os.path.exists(dirname):
+ print dirname
+# exit()
+ else:
+ os.mkdir(dirname)
+else:
+ p = 2
+ lam = 1e-8
+ use_mask = True
+ order = 1
+ dirname = "./"
-LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
logging.basicConfig(level=logging.INFO)
+
###################################################################################################
threshold = 1
a = 1.0 # in nm
@@ -38,11 +65,8 @@ inter = 'none'
dim = (1,) + (64, 64)
dim_small = (64, 64)
smoothed_pictures = True
-lam = 1E-4
-order = 1
log = True
PATH = './'
-dirname = PATH
###################################################################################################
# Read in files:
phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
@@ -50,16 +74,29 @@ phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
with open(PATH + 'mask.pickle') as pf:
mask = pickle.load(pf)
# Setup:
+if not use_mask:
+ mask = np.ones_like(mask, dtype=bool)
data_set = DataSet(a, dim, b_0, mask=mask)
data_set.append(phase_map, SimpleProjector(dim))
-regularisator = ZeroOrderRegularisator(lam)
-regularisator = FirstOrderRegularisator(mask, lam)
-print "OOO"
+if order == 0:
+ regularisator = ZeroOrderRegularisator(mask, lam, p)
+elif order == 1:
+ regularisator = FirstOrderRegularisator(mask, lam, p)
+else:
+ raise NotImplementedError
# Reconstruct the magnetic distribution:
tic = clock()
-mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator)
-#mag_data_rec = rc.optimize_nonlin(data_set, regularisator=regularisator)
+
+if p == 2:
+ mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator)
+else:
+ print regularisator.p
+ mag_data_rec = rc.optimize_nonlin(data_set, regularisator=regularisator)
+
+with open(dirname + "result.pickle", "wb") as pf:
+ import cPickle
+ cPickle.dump(mag_data_rec, pf)
# .optimize_simple_leastsq(phase_map, mask, b_0, lam=lam, order=order)
print 'reconstruction time:', clock() - tic
--
GitLab