Menu
▾
▴
[Pyphant-devel] SF.net SVN: pyphant:[672] trunk
|
From: <zk...@us...> - 2010-03-18 11:41:27
|
Revision: 672
http://pyphant.svn.sourceforge.net/pyphant/?rev=672&view=rev
Author: zklaus
Date: 2010-03-18 11:41:19 +0000 (Thu, 18 Mar 2010)
Log Message:
-----------
Merge branch 'master' into svn-trunk
* master:
Enh: Added maxima to mra anaylsis and incorporated it into the parameter estimation.
Modified Paths:
--------------
trunk/doc/demo/viewOSC.py
trunk/src/workers/OSC/OSC/ComputeFunctional.py
trunk/src/workers/OSC/OSC/EstimateParameter.py
trunk/src/workers/OSC/OSC/MRA.py
trunk/src/workers/OSC/OSC/tests/TestMRA.py
Modified: trunk/doc/demo/viewOSC.py
===================================================================
--- trunk/doc/demo/viewOSC.py 2010-03-14 14:26:07 UTC (rev 671)
+++ trunk/doc/demo/viewOSC.py 2010-03-18 11:41:19 UTC (rev 672)
@@ -171,7 +171,8 @@
worker = recipe.getWorker("ThicknessModeller")
simulation = worker.plugCalcAbsorption.getResult()
worker = recipe.getWorker("MRA Exp")
- minimaPos = worker.plugMra.getResult().inUnitsOf(simulation.dimensions[1])
+ minimaPos = worker.plugMra.getResult()[r'\lambda_{min}'].inUnitsOf(simulation.dimensions[1])
+ maximaPos = worker.plugMra.getResult()[r'\lambda_{max}'].inUnitsOf(simulation.dimensions[1])
worker = recipe.getWorker("AddColumn")
table = worker.plugCompute.getResult(subscriber=TextSubscriber("Add Column"))
xPos = table[u"x-position"]
@@ -196,6 +197,12 @@
label ="$\\Delta%s$"%minimaPos.shortname, linestyle='dashed')
pylab.vlines(minimaPos.data[:,index]-minimaPos.error[:,index],0.1,1.0,
label ="$\\Delta%s$"%minimaPos.shortname, linestyle='dashed')
+ pylab.vlines(maximaPos.data[:,index],0.1,1.0,
+ label ="$%s$"%maximaPos.shortname)
+ pylab.vlines(minimaPos.data[:,index]+maximaPos.error[:,index],0.1,1.0,
+ label ="$\\Delta%s$"%maximaPos.shortname, linestyle='dashed')
+ pylab.vlines(maximaPos.data[:,index]-maximaPos.error[:,index],0.1,1.0,
+ label ="$\\Delta%s$"%maximaPos.shortname, linestyle='dashed')
pylab.title(result)
pylab.xlabel(simulation.dimensions[1].label)
@@ -205,7 +212,8 @@
worker = recipe.getWorker("ThicknessModeller")[0]
simulation = worker.plugCalcAbsorption.getResult()
worker = recipe.getWorker("MRA Exp")
- minimaPos = worker.plugMra.getResult().inUnitsOf(simulation.dimensions[1])
+ minimaPos = worker.plugMra.getResult()[r'\lambda_{min}'].inUnitsOf(simulation.dimensions[1])
+ maximaPos = worker.plugMra.getResult()[r'\lambda_{max}'].inUnitsOf(simulation.dimensions[1])
worker = recipe.getWorker("AddColumn")
table = worker.plugCompute.getResult(subscriber=TextSubscriber("Add Column"))
xPos = table[u"x-position"]
@@ -222,6 +230,8 @@
if not noIndicators:
pylab.vlines(minimaPos.data[:,index],0.1,1.0,
label ="$%s$"%minimaPos.shortname)
+ pylab.vlines(maximaPos.data[:,index],0.1,1.0,
+ label ="$%s$"%maximaPos.shortname)
pylab.title(result)
pylab.xlabel(simulation.dimensions[1].label)
pylab.ylabel(simulation.label)
@@ -239,11 +249,12 @@
def functional(recipe, curvNo, noIndicators):
worker = recipe.getWorker("Compute Functional")
- functional = worker.plugCompute.getResult()
+ functional = worker.plugCompute.getResult()[r'F_{\lambda_{min}}']
worker = recipe.getWorker("ThicknessModeller")
simulation = worker.plugCalcAbsorption.getResult()
worker = recipe.getWorker("MRA Exp")
- minimaPos = worker.plugMra.getResult().inUnitsOf(simulation.dimensions[1])
+ minimaPos = worker.plugMra.getResult()[r'\lambda_{min}'].inUnitsOf(simulation.dimensions[1])
+ maximaPos = worker.plugMra.getResult()[r'\lambda_{max}'].inUnitsOf(simulation.dimensions[1])
worker = recipe.getWorker("AddColumn")
table = worker.plugCompute.getResult(subscriber=TextSubscriber("Add Column"))
thickness = table[u"thickness"]
@@ -253,6 +264,8 @@
ordinate = functional.dimensions[1].data
pylab.hlines(minimaPos.data[:,index],ordinate.min(),ordinate.max(),
label ="$%s$"%minimaPos.shortname)
+ pylab.hlines(maximaPos.data[:,index],ordinate.min(),ordinate.max(),
+ label ="$%s$"%maximaPos.shortname)
abscissae = functional.dimensions[0].data
pylab.vlines(thickness.data[index],abscissae.min(),abscissae.max(),
label ="$%s$"%minimaPos.shortname)
@@ -265,7 +278,8 @@
thickness = table[u"thickness"]
index = curvNo2Index(table[u"pixel"], curvNo)
worker = recipe.getWorker("MRA Exp")
- minimaPos = worker.plugMra.getResult().inUnitsOf(simulation.dimensions[1])
+ minimaPos = worker.plugMra.getResult()[r'\lambda_{min}'].inUnitsOf(simulation.dimensions[1])
+ maximaPos = worker.plugMra.getResult()[r'\lambda_{max}'].inUnitsOf(simulation.dimensions[1])
visualizer = ImageVisualizer(simulation, False)
ordinate = simulation.dimensions[1].data
if not noIndicators:
@@ -274,6 +288,8 @@
abscissae = simulation.dimensions[0].data
pylab.vlines(minimaPos.data[:,index],abscissae.min(),abscissae.max(),
label ="$%s$"%minimaPos.shortname)
+ pylab.vlines(maximaPos.data[:,index],abscissae.min(),abscissae.max(),
+ label ="$%s$"%maximaPos.shortname)
def dumpMinima(recipe):
worker = recipe.getWorker("ThicknessModeller")
@@ -286,7 +302,7 @@
thickness = table[u"thickness"]
cols = [index, xPos, yPos, thickness]
worker = recipe.getWorker("MRA Exp")
- minimaPos = worker.plugMra.getResult().inUnitsOf(simulation.dimensions[1])
+ minimaPos = worker.plugMra.getResult()[r'\lambda_{min}'].inUnitsOf(simulation.dimensions[1])
import numpy
data = numpy.vstack([ c.data for c in cols]
+ numpy.vsplit(minimaPos.data, len(minimaPos.data))).transpose()
Modified: trunk/src/workers/OSC/OSC/ComputeFunctional.py
===================================================================
--- trunk/src/workers/OSC/OSC/ComputeFunctional.py 2010-03-14 14:26:07 UTC (rev 671)
+++ trunk/src/workers/OSC/OSC/ComputeFunctional.py 2010-03-18 11:41:19 UTC (rev 672)
@@ -53,7 +53,7 @@
REVISION = "$Revision$"[11:-1]
name = "Compute Functional"
- _sockets = [("field", Connectors.TYPE_IMAGE)]
+ _sockets = [("field", Connectors.TYPE_ARRAY)]
_params = [("extentX", u"Extension of x-axis [%%]", 10, None),
("extentY", u"Extension of y-axis [%%]", 10, None)]
@@ -92,13 +92,15 @@
subscriber %= percentage
result = DataContainer.FieldContainer(functional.transpose(),
dimensions=[yDim, xDim],
- longname = 'functional',
- shortname= 'F'
+ longname = 'functional of %s'%field.longname,
+ shortname= 'F_{%s}'%field.shortname
)
return result
- @Worker.plug(Connectors.TYPE_IMAGE)
+ @Worker.plug(Connectors.TYPE_ARRAY)
def compute(self, field, subscriber=1):
- functional = self.computeDistances(field, subscriber)
- functional.seal()
- return functional
+ functionals = DataContainer.SampleContainer([self.computeDistances(column, subscriber) for column in field],
+ longname='Functionals of %s'%field.longname,
+ shortname='F_{%s}'%field.shortname)
+ functionals.seal()
+ return functionals
Modified: trunk/src/workers/OSC/OSC/EstimateParameter.py
===================================================================
--- trunk/src/workers/OSC/OSC/EstimateParameter.py 2010-03-14 14:26:07 UTC (rev 671)
+++ trunk/src/workers/OSC/OSC/EstimateParameter.py 2010-03-18 11:41:19 UTC (rev 672)
@@ -51,20 +51,34 @@
REVISION = "$Revision$"[11:-1]
name = "Estimate Parameter"
- _sockets = [("model", Connectors.TYPE_IMAGE),
- ("experimental", Connectors.TYPE_IMAGE)]
- _params = [("extentX", u"Extension of x-axis [%%]", 10, None),
+ _sockets = [("model", Connectors.TYPE_ARRAY),
+ ("experimental", Connectors.TYPE_ARRAY)]
+ _params = [("minima_model", u"Minima in the model", [u"Minima"], None),
+ ("maxima_model", u"Maxima in the model", [u"Maxima"], None),
+ ("minima_experimental", u"Minima in the experiment", [u"Minima"], None),
+ ("maxima_experimental", u"Maxima in the experiment", [u"Maxima"], None),
+ ("extentX", u"Extension of x-axis [%%]", 10, None),
("extentY", u"Extension of y-axis [%%]", 10, None)]
- def calculateThickness(self, row, model, error=None):
+ def refreshParams(self, subscriber=None):
+ if self.socketModel.isFull():
+ templ = self.socketModel.getResult( subscriber )
+ self.paramMinima_model.possibleValues = templ.longnames.keys()
+ self.paramMaxima_model.possibleValues = templ.longnames.keys()
+ if self.socketExperimental.isFull():
+ templ = self.socketExperimental.getResult( subscriber )
+ self.paramMinima_experimental.possibleValues = templ.longnames.keys()
+ self.paramMaxima_experimental.possibleValues = templ.longnames.keys()
+
+ def calculateThickness(self, row_minima, row_maxima, minima_model, maxima_model,
+ minima_error=None, maxima_error=None):
"""
Given a vector of minima (row) and a 2 dimensional model,
this estimates the corresponding parameter.
"""
- if len(row)==0:
+ if len(row_minima)+len(row_maxima)==0:
return numpy.nan
- data = model.data.transpose()
- def calc(row, col, error):
+ def calc(row, model, col, error):
if error:
weight=0
for c,e in zip(row, error):
@@ -78,38 +92,68 @@
else:
return sum([col[numpy.argmin((model.dimensions[0].data-c)**2)]
for c in row])
- i = numpy.argmin(numpy.array([calc(row, col, error) for col in data]))
- return model.dimensions[1].data[i]
+ minima_data = minima_model.data.transpose()
+ maxima_data = maxima_model.data.transpose()
+ i = numpy.argmin(numpy.array([calc(row_minima, minima_model, minima_col, minima_error)
+ +calc(row_maxima, maxima_model, maxima_col, maxima_error)
+ for minima_col, maxima_col in
+ zip(minima_data, maxima_data)]))
+ assert(minima_model.dimensions[1].data[i]==maxima_model.dimensions[1].data[i])
+ return minima_model.dimensions[1].data[i]
@Worker.plug(Connectors.TYPE_IMAGE)
def compute(self, model, experimental, subscriber=1):
- renormedExp = experimental.inUnitsOf(model.dimensions[0])
- minima = renormedExp.data.transpose()
- if renormedExp.error != None:
- error = iter(renormedExp.error.transpose())
+ minima_model = model[self.paramMinima_model.value]
+ maxima_model = model[self.paramMaxima_model.value]
+ minima_experimental = experimental[
+ self.paramMinima_experimental.value
+ ].inUnitsOf(
+ minima_model.dimensions[0]
+ )
+ maxima_experimental = experimental[
+ self.paramMaxima_experimental.value
+ ].inUnitsOf(
+ maxima_model.dimensions[0]
+ )
+ minima = minima_experimental.data.transpose()
+ if minima_experimental.error != None:
+ minima_error = iter(minima_experimental.error.transpose())
else:
- error = None
+ minima_error = None
+ maxima = maxima_experimental.data.transpose()
+ if maxima_experimental.error != None:
+ maxima_error = iter(maxima_experimental.error.transpose())
+ else:
+ maxima_error = None
parameter = []
inc = 100.0/float(len(minima))
acc = inc
subscriber %= acc
- for row in minima:
- if error:
- filteredError = filter(
- lambda c: not numpy.isnan(c), error.next())
+ for row_minima, row_maxima in zip(minima, maxima):
+ if minima_error:
+ filtered_minima_error = filter(
+ lambda c: not numpy.isnan(c), minima_error.next())
else:
- filteredError = None
+ filtered_minima_error = None
+ if maxima_error:
+ filtered_maxima_error = filter(
+ lambda c: not numpy.isnan(c), maxima_error.next())
+ else:
+ filtered_maxima_error = None
parameter.append(self.calculateThickness(
- filter(lambda c: not numpy.isnan(c), row),
- model,
- filteredError
+ filter(lambda c: not numpy.isnan(c), row_minima),
+ filter(lambda c: not numpy.isnan(c), row_maxima),
+ minima_model,
+ maxima_model,
+ filtered_minima_error,
+ filtered_maxima_error
))
acc += inc
subscriber %= acc
result = DataContainer.FieldContainer(
numpy.array(parameter),
- longname = model.dimensions[-1].longname,
- shortname = model.dimensions[-1].shortname,
- unit = model.dimensions[-1].unit)
+ longname = minima_model.dimensions[-1].longname,
+ shortname = minima_model.dimensions[-1].shortname,
+ unit = minima_model.dimensions[-1].unit)
result.seal()
return result
Modified: trunk/src/workers/OSC/OSC/MRA.py
===================================================================
--- trunk/src/workers/OSC/OSC/MRA.py 2010-03-14 14:26:07 UTC (rev 671)
+++ trunk/src/workers/OSC/OSC/MRA.py 2010-03-18 11:41:19 UTC (rev 672)
@@ -1,8 +1,9 @@
# -*- coding: utf-8 -*-
-# Copyright (c) 2008-2009, Rectorate of the University of Freiburg
+# Copyright (c) 2008-2010, Rectorate of the University of Freiburg
# Copyright (c) 2009, Andreas W. Liehr (li...@us...)
-# All rights reserved.
+# Copyright (c) 2010, Klaus Zimmermann (zk...@us...)
+# all rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
@@ -49,12 +50,29 @@
from scipy.special import ive
from scipy.signal import convolve
+def findMaxima(fieldData, lastExtrema=None):
+ leftLess = fieldData[:-2] < fieldData[1:-1]
+ leftEqual = fieldData[:-2] == fieldData[1:-1]
+ rightLess = fieldData[2:] < fieldData[1:-1]
+ rightEqual = fieldData[2:] == fieldData[1:-1]
+ maxima_c = numpy.logical_and(leftLess, rightLess)
+ maxima_le = numpy.logical_and(leftLess, rightEqual)
+ maxima_re = numpy.logical_and(rightLess, leftEqual)
+ maxima_e = numpy.logical_and(maxima_le[:-1], maxima_re[1:])
+ maxima = numpy.logical_or(maxima_c[1:], maxima_e).nonzero()[0]
+ if lastExtrema==None or len(maxima)==0 or len(lastExtrema)==len(maxima):
+ return maxima
+ trackedMaxima = []
+ for lastMaximum in lastExtrema:
+ distance = (maxima-lastMaximum)**2
+ trackedMaxima.append(distance.argmin())
+ return maxima[trackedMaxima]
+
def findMinima(fieldData, lastExtrema=None):
leftGreater = fieldData[:-2] > fieldData[1:-1]
leftEqual = fieldData[:-2] == fieldData[1:-1]
rightGreater = fieldData[2:] > fieldData[1:-1]
rightEqual = fieldData[2:] == fieldData[1:-1]
-
minima_c = numpy.logical_and(leftGreater, rightGreater)
minima_le = numpy.logical_and(leftGreater, rightEqual)
minima_re = numpy.logical_and(rightGreater, leftEqual)
@@ -75,17 +93,6 @@
kernel = ive(numpy.arange(-n, n), sigma)
return convolve(field.data, kernel, mode='same')
-#def mra1d(dim, field, n):
-# mrr = [convolveMRA(field, sigma) for sigma in
-# numpy.linspace(1, n,10).tolist()]
-# mrr.insert(0,field.data)
-# firstMinima = lastMinima = findMinima(mrr[-1], None)
-# for row in reversed(mrr[:-1]):
-# lastMinima = findMinima(row, lastMinima)
-# pos = dim.data[numpy.array(lastMinima)+1]
-# error = numpy.abs(pos - dim.data[numpy.array(firstMinima)+1])
-# return pos, error
-
class MraError(RuntimeError):
def __init__(self, msg, convolvedField):
RuntimeError.__init__(self, msg)
@@ -95,20 +102,30 @@
sigmaSpace = numpy.linspace(n, 1, 10)
convolvedField = convolveMRA(field, sigmaSpace[0])
firstMinima = lastMinima = findMinima(convolvedField, None)
- if len(firstMinima)==0:
- raise MraError("No Minima found at sigma level %s"%sigmaSpace[0],
+ firstMaxima = lastMaxima = findMaxima(convolvedField, None)
+ if len(firstMinima)==0 and len(firstMaxima)==0:
+ raise MraError("No Extrema found at sigma level %s"%sigmaSpace[0],
convolvedField)
for sigma in sigmaSpace[1:]:
convolvedField = convolveMRA(field, sigma)
lastMinima = findMinima(convolvedField, lastMinima)
- if len(lastMinima)==0:
- import pylab
- pylab.plot(convolvedField)
- pylab.show()
- pos = dim.data[numpy.array(lastMinima)+1]
- error = numpy.abs(pos - dim.data[numpy.array(firstMinima)+1])
- return pos, error
+ lastMaxima = findMaxima(convolvedField, lastMaxima)
+ pos_minima = dim.data[numpy.array(lastMinima)+1]
+ error_minima = numpy.abs(pos_minima - dim.data[numpy.array(firstMinima)+1])
+ pos_maxima = dim.data[numpy.array(lastMaxima)+1]
+ error_maxima = numpy.abs(pos_maxima - dim.data[numpy.array(firstMaxima)+1])
+ return ((pos_minima, error_minima), (pos_maxima, error_maxima))
+def pos_error_to_data_container(p_e):
+ n = max(map(lambda (p,e): len(p), p_e))
+ m = len(p_e)
+ pos = numpy.ones((m,n),'float')*numpy.NaN
+ error = pos.copy()
+ for i in xrange(m):
+ for j in xrange(len(p_e[i][0])):
+ pos[i,j] = p_e[i][0][j]
+ error[i,j] = p_e[i][1][j]
+ return n, pos, error
class MRA(Worker.Worker):
API = 2
@@ -121,7 +138,7 @@
("longname",u"Name of result",'default',None),
("symbol",u"Symbol of result",'default',None)]
- @Worker.plug(Connectors.TYPE_IMAGE)
+ @Worker.plug(Connectors.TYPE_ARRAY)
def mra(self, field, subscriber=0):
dim = field.dimensions[-1]
try:
@@ -139,30 +156,36 @@
try:
p_e.append(mra1d(dim, field1d, sigmaMax))
except MraError:
- p_e.append(([],[]))
+ p_e.append((([],[]),([],[])))
acc += inc
subscriber %= acc
- n = max(map(lambda (p,e): len(p), p_e))
- m = len(p_e)
- pos = numpy.ones((m,n),'float')*numpy.NaN
- error = pos.copy()
- for i in xrange(m):
- for j in xrange(len(p_e[i][0])):
- pos[i,j] = p_e[i][0][j]
- error[i,j] = p_e[i][1][j]
- dims = [DataContainer.generateIndex(0,n), field.dimensions[0]]
+ minima, maxima = zip(*p_e)
+ n_min, pos_min, err_min = pos_error_to_data_container(minima)
+ n_max, pos_max, err_max = pos_error_to_data_container(maxima)
+ dims_min = [DataContainer.generateIndex(0,n_min), field.dimensions[0]]
+ dims_max = [DataContainer.generateIndex(0,n_max), field.dimensions[0]]
else:
- pos, error = mra1d(dim, field, sigmaMax)
- n = len(pos)
- dims = [DataContainer.generateIndex(0,n)]
+ (pos_min, err_min), (pos_max, err_max) = mra1d(dim, field, sigmaMax)
+ dims_min = [DataContainer.generateIndex(0,len(pos_min))]
+ dims_max = [DataContainer.generateIndex(0,len(pos_max))]
subscriber %= 100.
- roots = DataContainer.FieldContainer(pos.transpose(),
- error = error.transpose(),
- unit = dim.unit,
- dimensions = dims,
- mask = numpy.isnan(pos).transpose(),
- longname="%s of the local %s of %s" % (dim.longname,"minima",field.longname),
- shortname="%s_0" % dim.shortname)
+ minima = DataContainer.FieldContainer(pos_min.transpose(),
+ error = err_min.transpose(),
+ unit = dim.unit,
+ dimensions = dims_min,
+ mask = numpy.isnan(pos_min).transpose(),
+ longname="%s of the local %s of %s" % (dim.longname,"minima",field.longname),
+ shortname="%s_{min}" % dim.shortname)
+ maxima = DataContainer.FieldContainer(pos_max.transpose(),
+ error = err_max.transpose(),
+ unit = dim.unit,
+ dimensions = dims_max,
+ mask = numpy.isnan(pos_max).transpose(),
+ longname="%s of the local %s of %s" % (dim.longname,"maxima",field.longname),
+ shortname="%s_{max}" % dim.shortname)
+ roots = DataContainer.SampleContainer([minima, maxima],
+ longname="%s of the local %s of %s" % (dim.longname,"extrema",field.longname),
+ shortname="%s_{extrem}" % dim.shortname)
if self.paramLongname.value != 'default':
roots.longname = self.paramLongname.value
if self.paramSymbol.value != 'default':
Modified: trunk/src/workers/OSC/OSC/tests/TestMRA.py
===================================================================
--- trunk/src/workers/OSC/OSC/tests/TestMRA.py 2010-03-14 14:26:07 UTC (rev 671)
+++ trunk/src/workers/OSC/OSC/tests/TestMRA.py 2010-03-18 11:41:19 UTC (rev 672)
@@ -87,8 +87,26 @@
w.paramScale.value = "1.0m"
result = w.mra(self.V)
#Testing
- numpy.testing.assert_array_almost_equal(result.data,expectedResult.data,4)
+ numpy.testing.assert_array_almost_equal(result[r'x_{min}'].data,expectedResult.data,4)
+ def testMaxima(self):
+ """Test the correct computation of all local minima for a bistable potential."""
+ #Predict result
+ x0,curv,mask = fixedPoints(numpy.array([self.LAMBDA]),kappa1=self.kappa1)
+ expectedResult = DC.FieldContainer(numpy.extract(curv[0]>0,x0[0]),
+ unit = self.xField.unit,
+ longname = 'position of the local minima of electric potential',
+ shortname = 'x_0')
+ #Retrieve result from worker
+ w = MRA.MRA(None)
+ w.paramScale.value = "1.0m"
+ V = copy.deepcopy(self.V)
+ V.data*=-1
+ result = w.mra(V)
+ #Testing
+ numpy.testing.assert_array_almost_equal(result[r'x_{max}'].data,expectedResult.data,4)
+
+
class TestExtremumFinderTable(unittest.TestCase):
"""Sets up a mirror symmetric bistable potential with a continuous
distretisation and computes its local extrema and the respective
@@ -142,7 +160,7 @@
w = MRA.MRA(None)
w.paramScale.value = "1.0m"
#Retrieve result from worker
- result = copy.deepcopy(w.mra(self.V))
+ result = copy.deepcopy(w.mra(self.V))['x_{min}']
result.error=None
self.test(result,expectedResult,1e-2,1e-2)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
