Archief - [Python] The truth value of an array...

Exorikos · 21 okt 2012

Ik zou graag een chisqr-plot maken.

Eerst fit ik zonder problemen een spectrum aan de experimentele waarden. Daarna wil ik de twee belangrijke parameters die hieruit komen varieren en voor elke variatie de chisqr berekenen. Hij geeft echter een error met traceback naar een stuk code dat daarvoor volledig werkte: ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Hieronder vind je alle relevante stukken uit de code. De traceback gaat naar "if E < centroid:", een boolean die niet ambigu is en waar ik niet onderuit kan.

Code:

# Builds spectrum
def spectrum(x,al,centroid,bkg,I,J,HWHM,intensityl,intensityr):
	F = states(I,J)
	build = zeros(len(x))
	for i in F:
		E = centroid - energy(i,I,J,al)
		if E < centroid:
			build+=lorentzian(x,E,HWHM,intensityl)
		else:
			build+=lorentzian(x,E,HWHM,intensityr)
	build+=bkg
	return build

# Function that calculates the chi square for the chi square map

def chisquare(x,y,al,centroid,bkg,I,J,hwhm,intensityl,intensityr):
	y = array(y)
	expected = spectrum(x,al,centroid,bkg,I,J,hwhm,intensityl,intensityr)
	chisq = 0
	for i in range(len(y)):
		chisq+=(y[i]-expected[i])**2/expected[i]
	return chisq

# Script voor het chisqr plot
com = linspace(params['centroid'].value-1000,params['centroid'].value+1000,20)
a = linspace(params['al'].value-500,params['centroid'].value+500,10)

COM,A = meshgrid(com,a)
contourf(COM,A,chisquare(x,y,A,COM,params['bkg'].value,I,J,params['HWHM'].value,params['intensityl'].value,params['intensityr'].value))
show()

Hopelijk is het wat duidelijk. Extra info verstrek ik graag.

Fraggie · 21 okt 2012

Daar Python aan Duck Typing doet valt er weinig te zeggen over dat stukje code tenzij je aangeeft wat energy(..) terug geeft en wat een centroid is.

Afhankelijk van "hoe" object georiënteerd je werkt kan wel wat aan type hinting (is geen feature van de taal!) doen, door het type wél expliciet te specificeren.

Een andere oplossing kan zijn, om na de bewerking centroid - energy(...) het resultaat expliciet te converteren naar een type eigen aan Python (int, long, float, complex..) of het type van de library dat je gebruikt (zoals bv numpy/scipy).

Exorikos · 21 okt 2012

Normaal zouden al die dingen floats moeten zijn. Ik heb ze nu als float gecast. Omdat ik in die contourf al de waarden als array doorgeef kan Python er blijkbaar niet mee om. Mijn plotcommando is dus fout? Is er een andere manier?

Fraggie · 21 okt 2012

Hangt er vanaf, daarvoor moet ik meer code zien, maar neem bv het volgende simpele voorbeeld met NumPy:

Code:

>>> a = np.arange(12).reshape(3,4)
>>> a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> a > 2
array([[False, False, False,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True]], dtype=bool)
>>> if (a > 2): print "troll"
else: print "notroll"


Traceback (most recent call last):
  File "<pyshell#64>", line 1, in <module>
    if (a > 2): print "troll"
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

de operator > heeft een ander effect afhankelijk van de datastructuur die gebruikt wordt.

Messias. · 21 okt 2012

Bon ja, die operator geeft dus voor die twee datatypes een niet-scalair resultaat terug. Kunt ge natuurlijk nie gebruiken als input voor een if-statement. Wat moet dat if-statement aanvangen met een array van true's en falsen?

Vandaar dat de foutboodschap ook suggereert om any() of all() te gebruiken om de collectie van booleans te aggregeren in één waarde.

Exorikos · 21 okt 2012

Code:

from numpy import *
from pylab import *
from scipy import constants
import glob
from lmfit import minimize, Parameters, Parameter, report_errors
from operator import *
import os.path
from Tkinter import *
from mpl_toolkits.mplot3d import Axes3D

# Define variables and parameters

I = 4.5
J = 0.5
iscool = 49970.0
mass = 221.014254

params = Parameters()
params.add('al', value=4000.0)
params.add('centroid', value=3000.0)
params.add('HWHM', value=800.0,min=0.0)
params.add('bkg', value=0.0)
params.add('intensityl', value=700.0, min=0.0)
params.add('intensityr', value=400.0, min=0.0)

# Working out the different couplings of I and J

def states(I,J):
	F = list()
	value = abs(I-J)
	while value < I+J+1:
		F.append(value)
		value += 1
	F = array(F)
	return F

# Calculating C

def calculateC(F,I,J):
	C = F*(F+1) - I*(I+1) - J*(J+1)
	return C

# Calculating the energy difference to the unperturbed state

def energy(F,I,J,al):
	C = calculateC(F,I,J)
	E = -C*al/2
	return E

# Builds a Lorentzian peak

def lorentzian(x,x0,HWHM,intensity):
	a = intensity * (HWHM**2/((x-x0)**2+HWHM**2)) #Lorentz-Cauchy distr Wikipedia
	return a

# Builds spectrum

def spectrum(x,al,centroid,bkg,I,J,HWHM,intensityl,intensityr):
	F = states(I,J)
	build = zeros(len(x))
	for i in F:
		E = float(centroid - energy(i,I,J,al))
		if E < centroid:
			build+=lorentzian(x,E,HWHM,intensityl)
		else:
			build+=lorentzian(x,E,HWHM,intensityr)
	build+=bkg
	return build

# Doppler shift through ISCOOL

def doppler(frequency):
	beta = sqrt(1 - (mass**2*constants.c**4)/(iscool+mass*constants.c**2)**2)
	freq = frequency*(1-beta)/sqrt(1-beta**2)
	return freq

# Convert wavenumber to frequency

def convert(wavenumber):
	freq = 2*100*constants.c*wavenumber
	freq = doppler(freq)
	freq = freq/1000000
	freq = freq - 709726866.314
	return freq

# Defines the difference between the expected and the experimental spectrum

def residuals(params,x,y):
	al = params['al'].value
	centroid = params['centroid'].value
	hwhm = params['HWHM'].value
	intensityl = params['intensityl'].value
	intensityr = params['intensityr'].value
	bkg = params['bkg'].value
	y = array(y)
	expected = spectrum(x,al,centroid,bkg,I,J,hwhm,intensityl,intensityr)
	err = map(sub,y,expected)
	return err/sqrt(y+1)

# Function that calculates the chi square for the chi square map

def chisquare(x,y,al,centroid,bkg,I,J,hwhm,intensityl,intensityr):
	y = array(y)
	expected = spectrum(x,al,centroid,bkg,I,J,hwhm,intensityl,intensityr)
	chisq = 0
	for i in range(len(y)):
		chisq+=(y[i]-expected[i])**2/expected[i]
	return chisq


# Import frequency data file from directory

for filename in glob.glob("*.txt"):


# Open file and import frequencies and counts into arrays

	data = open(filename, "r")
	lines = data.readlines()
	x=zeros(len(lines))
	wavenumber=zeros(len(lines))
	y=zeros(len(lines))
	n=0
	for line in lines:
			newline = line.split()
			xdat=newline[0]
			ydat=newline[1]
			x[n]=float(xdat);y[n]=float(ydat)
			wavenumber[n]=x[n]
			x[n]=convert(x[n])	
			n=n+1




# Do fit
	results = minimize(residuals, params, args=(x,y),maxfev=1000000000)

# 
	com = linspace(params['centroid'].value-1000,params['centroid'].value+1000,20)
	a = linspace(params['al'].value-500,params['centroid'].value+500,10)

	COM,A = meshgrid(com,a)
	contourf(COM,A,chisquare(x,y,A,COM,params['bkg'].value,I,J,params['HWHM'].value,params['intensityl'].value,params['intensityr'].value))
	show()

De (bijna) volledige code. De error die hij geeft nu ik cast naar float is dat hij geen array groter dan lengte 1 kan converteren naar float. Kortom, Python kan niet overweg met de arrays die ik doorstuur vanuit dat contourf commando. Ik had verwacht dat Python ermee overweg kon, want ik had een voorbeeld gevonden waarin stond: contourf(X,Y,f(X,Y))

EDIT: Hopelijk is mijn code niet vreselijk ongestructureerd...

EDIT2: Het lijkt alsof hij bij energy in plaats van een waarde voor "al" hij ineens het hele array van de variaties doorgeeft. Daardoor geeft hij daar een array terug in plaats van een float.

iterums · 22 okt 2012

In contourf gebruik je als derde argument

Code:

chisquare(x,y,A,COM,params['bkg'].value,I,J,params['HWHM'].value,params['intensityl'].value

nadat

Code:

COM,A = meshgrid(com,a)

Je gebruikt dus als derde argument 'al' voor chisquare 'A', wat uiteraard niet in overeenstemming is met wat je verwacht, ook al ik heb geen idee van wat dat nu precies is (-1 voor onduidelijke variabelenamen

).

Zie documentatie meshgrid:

Code:

def meshgrid(x,y):
    """
    Return coordinate matrices from two coordinate vectors.

    Parameters
    ----------
    x, y : ndarray
        Two 1-D arrays representing the x and y coordinates of a grid.

    Returns
    -------
    X, Y : ndarray
        For vectors `x`, `y` with lengths ``Nx=len(x)`` and ``Ny=len(y)``,
        return `X`, `Y` where `X` and `Y` are ``(Ny, Nx)`` shaped arrays
        with the elements of `x` and y repeated to fill the matrix along
        the first dimension for `x`, the second for `y`.
(...)

Duidelijk definiëren van verwachte invoerargumenten en uitvoer...

Exorikos · 22 okt 2012

iterums zei:
In contourf gebruik je als derde argument

Code:

chisquare(x,y,A,COM,params['bkg'].value,I,J,params['HWHM'].value,params['intensityl'].value

nadat

Code:

COM,A = meshgrid(com,a)

Je gebruikt dus als derde argument 'al' voor chisquare 'A', wat uiteraard niet in overeenstemming is met wat je verwacht, ook al ik heb geen idee van wat dat nu precies is (-1 voor onduidelijke variabelenamen ).

Zie documentatie meshgrid:

Code:

def meshgrid(x,y): """ Return coordinate matrices from two coordinate vectors. Parameters ---------- x, y : ndarray Two 1-D arrays representing the x and y coordinates of a grid. Returns ------- X, Y : ndarray For vectors `x`, `y` with lengths ``Nx=len(x)`` and ``Ny=len(y)``, return `X`, `Y` where `X` and `Y` are ``(Ny, Nx)`` shaped arrays with the elements of `x` and y repeated to fill the matrix along the first dimension for `x`, the second for `y`. (...)

Duidelijk definiëren van verwachte invoerargumenten en uitvoer...

Is het niet meteen duidelijk dat die a variabele staat voor de hyperfijnsplitting door het magnetisch moment van de atoomkern?

Ik programmeer dit natuurlijk voor mensen in het veld, dat maakt het voor jullie inderdaad lastig.

Het is dus de bedoeling van de twee meest relevante parameters (center of mass en a) uit een model te varieren en te plotten. Alles wat die chisquare methode nodig heeft zijn de experimentele punten (x,y) en de modelparameters (com, a, bkg, hwhm, I, J, intensityl, intensityr). Nu heb ik twee for loops ingevoegd in chisquare:

Code:

# Function that calculates the chi square for the chi square map

def chisquare(x,y,al,centroid,bkg,I,J,hwhm,intensityl,intensityr):
	y = array(y)
	chisq=zeros((len(centroid),len(al)))
	l=0
	for k in centroid:
		m=0
		for n in al:
			expected = spectrum(x,n,k,bkg,I,J,hwhm,intensityl,intensityr)
			chisqr=0
			for i in range(len(y)):
				chisqr+=(y[i]-expected[i])**2/expected[i]
			chisq[l][m]=chisqr
			m+=1
		l+=1
	return chisq
# Script for naking chi**2 plot

	com = linspace(params['centroid'].value-1000,params['centroid'].value+1000,50)
	a = linspace(params['al'].value-500,params['centroid'].value+500,20)
	chisquared = chisquare(x,y,a,com,params['bkg'].value,I,J,params['HWHM'].value,params['intensityl'].value,params['intensityr'].value)


	COM,A = meshgrid(com,a)
	contourf(COM,A,chisquared)
	show()

Zo krijg ik toch een 2D arrray terug dacht ik, maar contourf klaagt dat hij een 1D of 2D array verwacht voor x en y. Dat krijgt hij toch?

Even geprint en hij krijgt een 2D array voor X,Y en een 2D array voor Z...

EDIT: Vreemd dat het wel werkt als ik x en y even grote arrays geef... Ik had nergens zien staan dat dat een voorwaarde was.

EDIT2: Nu nog uitzoeken hoe ik Python het script parallel voor meerdere datafiles laat lopen over 3 van de 4 cores.

Fraggie · 22 okt 2012

Exorikos zei:
EDIT2: Nu nog uitzoeken hoe ik Python het script parallel voor meerdere datafiles laat lopen over 3 van de 4 cores.

Zeker dat dit nodig is? Het array type van NumPy brengt het zogenaamde SIMD taxonomy in spel. Dus zonder dat je iets speciaal doet gebruikt deze maximale resources. Heb je 4 cores, dan gebruikt ie er 4, heb je er 12, dan gebruikt ie er standaard 12.

Wanneer je dit manueel gaat willen doen (en veel gaat vloeken op de GIL) kan je veel performance verliezen.

Vergeet ook niet dat het inlezen van files m.b.v. NumPy véél sneller is dan het inlezen met een andere module.

Exorikos · 22 okt 2012

Het is gewoon iets extra. Het is niet nodig, maar nu gebruikt hij maar één core voor het berekenen van de chisq map (berekening hiervan duurt het langst). Het zou leuk zijn als je elke datafile dus een aparte core kon geven. Drie (nog een core overlaten voor gewoon gebruik) files tegelijk laten analyseren is gewoon handig als je een hondertal files moet doen.

Fraggie · 23 okt 2012

Als hij maar één core gebruikt, dan zit er nog te veel klassieke Python code in je bewerkingen. Wanneer je NumPy begint te gebruiken heb je bijna geen enkele for meer nodig. Vergeet ook zo veel mogelijk die len(), klassieke loops, gebruik xrange() over range() in geneste for's.

Probeer zo veel mogelijk met matrices en arrays te werken. De "overhead" van het opstellen van deze matrices/arrays wordt te niet gedaan door de snelheid van NumPy. In mijn eindwerk van +6K lijnen worden 14 fors gebruikt, waarvan 8 in GUI, 2 voor design patterns, 2 voor een generic plotter in matplotlib en 2 in code.

Werk je echter zelf multithreaded (want multicore programmeren is nog iets anders, en ik weet niet of dit mogelijk is met Python + NumPy), dan ga je te veel overhead hebben o.w.v. de message passing etc. Zoals ik eerder al zei, in Python kan je standaard niet o.w.v. de GIL parallel schrijven naar één datastructuur vanuit verschillende threads. Je moet dus manueel al je datastructuren opdelen, afhankelijk van hoeveel threads je wil, verwerken (zonder garantie op parallellisme, alleen concurrent behaviour staat in code) en ze manueel mergen. Het is dan aan jou OS om te beslissen wat er parallel zal gebeuren.

blackrabbit · 23 okt 2012

Fraggie zei:
multithreaded (want multicore programmeren is nog iets anders,..),

Dat mag je toch effe verduidelijken

Fraggie · 23 okt 2012

blackrabbit zei:
Dat mag je toch effe verduidelijken

Klassiek multithreadded programmeren geeft de mogelijkheid om parallel uitgevoerd te worden, maar dat moet niet (= concurrency). 1 CPU kan toch 50 threads hebben.

Echt parallel programmeren, zoals ik gedaan heb in OpenCL op mijn GPU en CPU, spreek niet zo zeer over threads, maar introduceert tal van nieuwe jargon. Zo kan je met OpenCL kiezen om bv, parallel op 128 (simpele) cores te rekenen in zogenaamde kernels (+- het kleinste deeltje van je algoritme dat opgelost kan worden). Dit is wat men verstaat onder parllelisme en ikzelf in de eerder uitleg doortrek naar multicore programmeren.

Maar ik ben zeker dat je dat al wist

.

blackrabbit · 24 okt 2012

Mja, je gebruikte gewoon rare termen, denk ik.

Er is gewoon een verschil tussen concurrency en paralellism. Los van CPU/GPU/multicore/manycore. Parallelism draait om opvoeren van execution speed (je doet meer werk in dezelfde tijd), terwijl concurrency eerder draait om responsiveness (vb: UI) of inherente gelijktijdige operaties (vb: bank-transacties).

GPU-programming is inderdaad enkel nuttig voor parallellisme, maar het is GPU-programming dat dat nieuwe jargon introduceert.

Ook: vanuit het standpunt van de programmeur maakt het weinig uit of een multi-threaded programma op een single- dan wel multi-core gaat draaien. Zoals je zegt gaat het systeem (VM/OS/..) de verschillende threads tijd geven en zit je met dezelfde problemen (vb racing conditions) in beide gevallen.

Fraggie · 24 okt 2012

Die termen wil ik nu nog in de midden laten.

Ons geschil ligt hem waarschijnlijk vanuit welk oogpunt we er naar kijken. Jou uitleg lijkt me dan ook meer een eind resultaat van concurrency en parallelism. Daar waar ik kijk uit oogpunt van de hardware & software.

Stel nu dat ik een FPGA ontwerp in VHDL dan ontwerp ik mijn "paralelle hardware" en codeer ik mijn interacties concurrent. Dit staat dan bv los van de term responsiveness (digitale logica is gewoon digitaal, het werkt of het werkt niet). Ook moet mijn execution speed er niet op vooruit gaan als ik nu 32 flipflops implementeer of er maar 10. Een FPGA op 1MHz blijft op 1MHz werken (beetje zoals een signaal van 100Mbps even snel propageert als een van 1Gbps).

Archief - [Python] The truth value of an array...

Exorikos

Legacy Member

Fraggie

Legacy Member

Exorikos

Legacy Member

Fraggie

Legacy Member

Messias.

Legacy Member

Exorikos

Legacy Member

iterums

Legacy Member

Exorikos

Legacy Member

Fraggie

Legacy Member

Exorikos

Legacy Member

Fraggie

Legacy Member

blackrabbit

Legacy Member

Fraggie

Legacy Member

blackrabbit

Legacy Member

Fraggie

Legacy Member