Mercurial > lcfOS
changeset 97:5a965e9664f2
Movage
| author | windel |
|---|---|
| date | Mon, 24 Dec 2012 13:32:54 +0100 |
| parents | a350055d6119 |
| children | 3f772feb12ef |
| files | python/__init__.py python/astviewer.py python/bouncing_cube_opencl.py python/codeeditor.py python/integration_methods.py python/libs/project.py python/libs/widgets/__init__.py python/libs/widgets/astviewer.py python/libs/widgets/codeeditor.py python/project.py |
| diffstat | 10 files changed, 826 insertions(+), 215 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/__init__.py Mon Dec 24 13:32:54 2012 +0100 @@ -0,0 +1,5 @@ +# File to make this directory a package. + +from .codeeditor import CodeEdit +from .astviewer import AstViewer +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/astviewer.py Mon Dec 24 13:32:54 2012 +0100 @@ -0,0 +1,36 @@ +from PyQt4.QtCore import * +from PyQt4.QtGui import * + +def astToNamedElement(astNode, parentNode): + """ Helper to convert and AST tree to NamedElement tree: """ + item = QStandardItem(str(astNode)) + item.setData(astNode) + parentNode.appendRow(item) + for c in astNode.getChildren(): + astToNamedElement(c, item) + +# The actual widget: +class AstViewer(QTreeView): + sigNodeSelected = pyqtSignal(object) + def __init__(self, parent=None): + super(AstViewer, self).__init__(parent) + self.setHeaderHidden(True) + self.clicked.connect(self.selectHandler) + + def setAst(self, ast): + """ Create a new model and add all ast elements to it """ + model = QStandardItemModel() + if ast: + astToNamedElement(ast, model.invisibleRootItem()) + self.setModel( model ) + self.expandAll() + + def selectHandler(self, index): + if not index.isValid(): + return + model = self.model() + item = model.itemFromIndex(index) + node = item.data() + self.sigNodeSelected.emit(node) + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/bouncing_cube_opencl.py Mon Dec 24 13:32:54 2012 +0100 @@ -0,0 +1,406 @@ +from PyQt4.QtGui import * +from PyQt4.QtCore import * +from PyQt4.QtOpenGL import QGLWidget +from OpenGL.GL import * +from OpenGL.GLU import gluPerspective +import sys +from random import random +from math import pi, cos, sin, fabs, sqrt +from numpy import mat, array, ones, zeros, eye +from numpy.linalg import norm +import numpy as np +import time +import scipy.integrate +import pyopencl + +""" + Test script that lets a dice bounce. + Converted from 20-sim equations into python code. + + 20-sim website: + http://www.20sim.com + +""" +def drawCube(w): + glBegin(GL_QUADS) # Start Drawing The Cube + glColor3f(0.0,1.0,0.0) # Set The Color To Blue + glVertex3f( w, w,-w) # Top Right Of The Quad (Top) + glVertex3f(-w, w,-w) # Top Left Of The Quad (Top) + glVertex3f(-w, w, w) # Bottom Left Of The Quad (Top) + glVertex3f( w, w, w) # Bottom Right Of The Quad (Top) + + glColor3f(1.0,0.5,0.0) # Set The Color To Orange + glVertex3f( w,-w, w) # Top Right Of The Quad (Bottom) + glVertex3f(-w,-w, w) # Top Left Of The Quad (Bottom) + glVertex3f(-w,-w,-w) # Bottom Left Of The Quad (Bottom) + glVertex3f( w,-w,-w) # Bottom Right Of The Quad (Bottom) + + glColor3f(1.0,0.0,0.0) # Set The Color To Red + glVertex3f( w, w, w) # Top Right Of The Quad (Front) + glVertex3f(-w, w, w) # Top Left Of The Quad (Front) + glVertex3f(-w,-w, w) # Bottom Left Of The Quad (Front) + glVertex3f( w,-w, w) # Bottom Right Of The Quad (Front) + + glColor3f(1.0,1.0,0.0) # Set The Color To Yellow + glVertex3f( w,-w,-w) # Bottom Left Of The Quad (Back) + glVertex3f(-w,-w,-w) # Bottom Right Of The Quad (Back) + glVertex3f(-w, w,-w) # Top Right Of The Quad (Back) + glVertex3f( w, w,-w) # Top Left Of The Quad (Back) + + glColor3f(0.0,0.0,1.0) # Set The Color To Blue + glVertex3f(-w, w, w) # Top Right Of The Quad (Left) + glVertex3f(-w, w,-w) # Top Left Of The Quad (Left) + glVertex3f(-w,-w,-w) # Bottom Left Of The Quad (Left) + glVertex3f(-w,-w, w) # Bottom Right Of The Quad (Left) + + glColor3f(1.0,0.0,1.0) # Set The Color To Violet + glVertex3f( w, w,-w) # Top Right Of The Quad (Right) + glVertex3f( w, w, w) # Top Left Of The Quad (Right) + glVertex3f( w,-w, w) # Bottom Left Of The Quad (Right) + glVertex3f( w,-w,-w) # Bottom Right Of The Quad (Right) + glEnd() # Done Drawing The Quad + +def drawFloor(w, h): + glBegin(GL_QUADS) # Start Drawing The Cube + + glColor3f(1.0,0.5,0.0) # Set The Color To Orange + glVertex3f( w,-w,h)# Top Right Of The Quad (Bottom) + glVertex3f(-w,-w,h)# Top Left Of The Quad (Bottom) + glVertex3f(-w,w,h)# Bottom Left Of The Quad (Bottom) + glVertex3f( w,w,h)# Bottom Right Of The Quad (Bottom) + glEnd() # Done Drawing The Quad + +def drawAxis(): + glLineWidth(0.5) + glBegin(GL_LINES) + glColor3f(1.0, 0.0, 0.0) + glVertex3f(0,0,0) + glVertex3f(1,0,0) + glColor3f(0.0, 1.0, 0.0) + glVertex3f(0,0,0) + glVertex3f(0,1,0) + glColor3f(0.0, 0.0, 1.0) + glVertex3f(0,0,0) + glVertex3f(0,0,1) + glEnd() + + +def cross(A, B): + a = A.A1 + b = B.A1 + return mat(np.cross(a, b)).T + +def skew(X): + Y = mat(zeros( (3, 3) )) + a,b,c = X.A1 + Y[0,1] = -c + Y[0,2] = b + Y[1,0] = c + Y[1,2] = -a + Y[2,0] = -b + Y[2,1] = a + return Y + +def adjoint(T): + W = T[0:3, 0] + V = T[3:6, 0] + a = mat(zeros( (6,6) ) ) + a[0:3, 0:3] = skew(W) + a[3:6, 0:3] = skew(V) + a[3:6, 3:6] = skew(W) + return a + +def Adjoint(H): + R = H[0:3, 0:3] + P = H[0:3, 3] + a = mat(zeros( (6,6) ) ) + a[0:3, 0:3] = R + a[3:6, 3:6] = R + a[3:6, 0:3] = skew(P) * R + return a + +def quatToR(q): + x, y, z, w = q.A1 + r = mat(eye(3)) + r[0,0] = 1 - (2*y**2+2*z**2) + r[0,1] = 2*x*y+2*z*w + r[0,2] = 2*x*z - 2*y*w + r[1,0] = 2*x*y-2*z*w + r[1,1] = 1 - (2*x**2 + 2*z**2) + r[1,2] = 2*y*z + 2*x*w + r[2,0] = 2*x*z+2*y*w + r[2,1] = 2*y*z - 2*x*w + r[2,2] = 1 - (2*x**2+2*y**2) + return r + +def rotateAbout(axis, angle): + ax, ay, az = (axis/norm(axis)).A1 + qx = ax*sin(angle/2.0) + qy = ay*sin(angle/2.0) + qz = az*sin(angle/2.0) + qw = cos(angle/2.0) + q = mat(array([qx,qy,qz,qw])).T + return q + +def normalizeQuaternion(quat): + x,y,z,w = quat.A1 + magnitude = sqrt(x*x + y*y + z*z + w*w) + x = x / magnitude + y = y / magnitude + z = z / magnitude + w = w / magnitude + quat[0, 0] = x + quat[1, 0] = y + quat[2, 0] = z + quat[3, 0] = w + return quat + +def VTo4x4(V): + v1, v2, v3 = V.A1 + return mat(array( \ + [[0.0, -v3, v2, -v1], \ + [ v3, 0.0, -v1, -v2], \ + [-v2, v1, 0.0, -v3], \ + [v1, v2, v3, 0.0] ])) + +def homogeneous(R,p): + H = mat(eye(4)) + H[0:3, 0:3] = R + H[0:3, 3] = p + return H + +def rateOfChange(states, thetime, parameters): + quat = states[0:4, 0] # Orientation (4) + pos = states[4:7, 0] # Position (3) + P = states[7:13, 0] # Momentum (6) + massI, gravity = parameters + # Rigid body parts: + # Forward Kinematic chain: + H = homogeneous(quatToR(quat), pos) # Forward kinematics + + AdjX2 = mat(eye(6)) # The connectionpoint in the real world + adjAdjX2_1 = adjoint(AdjX2[0:6,0]) + adjAdjX2_2 = adjoint(AdjX2[0:6,1]) + adjAdjX2_3 = adjoint(AdjX2[0:6,2]) + adjAdjX2_4 = adjoint(AdjX2[0:6,3]) + adjAdjX2_5 = adjoint(AdjX2[0:6,4]) + adjAdjX2_6 = adjoint(AdjX2[0:6,5]) + AdjInv2 = Adjoint(H.I) + M2 = AdjInv2.T * (massI * AdjInv2) # Transfor mass to base frame + MassMatrix = M2 + + wrenchGrav2 = mat( zeros((1,6)) ) + wrenchGrav2[0, 0:3] = -cross(gravity, pos).T + wrenchGrav2[0, 3:6] = gravity.T + + Bk = mat( zeros( (6,6) )) + Bk[0:3, 0:3] = skew(P[0:3, 0]) + Bk[3:6, 0:3] = skew(P[3:6, 0]) + Bk[0:3, 3:6] = skew(P[3:6, 0]) + + # TODO: do this a cholesky: + v = np.linalg.solve(MassMatrix, P) # Matrix inverse like thingy ! + + T2_00 = v # Calculate the relative twist! + TM2 = T2_00.T * M2 + twistExternal = T2_00 # Twist van het blokje + TMSum2 = TM2 + + PDotBodies = mat( zeros( (6,1)) ) + PDotBodies[0,0] = TMSum2 * (adjAdjX2_1 * T2_00) + PDotBodies[1,0] = TMSum2 * (adjAdjX2_2 * T2_00) + PDotBodies[2,0] = TMSum2 * (adjAdjX2_3 * T2_00) + PDotBodies[3,0] = TMSum2 * (adjAdjX2_4 * T2_00) + PDotBodies[4,0] = TMSum2 * (adjAdjX2_5 * T2_00) + PDotBodies[5,0] = TMSum2 * (adjAdjX2_6 * T2_00) + + PDot = -PDotBodies - Bk * v + PDot += wrenchGrav2.T + + ##### Contact wrench part: + HB_W = H # Is H-matrix van het blokje + WrenchB = mat(zeros( (1,6) )) + for px in [-0.5, 0.5]: + for py in [-0.5, 0.5]: + for pz in [-0.5, 0.5]: + HB1_B = homogeneous(mat(eye(3)), mat([px,py,pz]).T) + HB1_W = HB_W * HB1_B + HW1_W = homogeneous(mat(eye(3)), HB1_W[0:3,3]) + HW_W1 = HW1_W.I + HB_W1 = HW_W1 * HB_W + + AdjHB_W1 = Adjoint(HB_W1) + TB_W1_W1 = AdjHB_W1 * twistExternal + z = HB1_W[2, 3] + vx, vy, vz = TB_W1_W1[3:6, 0].A1 + if z < 0: + # Contact forces: + Fx = -50.0*vx + Fy = -50.0*vy + Fz = -z*50000.0 + else: + Fx = 0.0 + Fy = 0.0 + Fz = 0.0 + # TODO: reflect impulse + WrenchW1 = mat([0,0,0,0,0,Fz]) + # Transform it back: + WrenchB += (AdjHB_W1.T * WrenchW1.T).T + ##### End of contact wrench + + PDot += (WrenchB * AdjInv2).T + + # Position and orientation rates: + QOmega = VTo4x4(v[0:3, 0]) + quatDot = 0.5 * QOmega * quat + vel = v[3:6, 0] + posDot = skew(v[0:3]) * pos + vel + # The rate vector: + rates = mat(zeros( (13,1) )) + rates[0:4, 0] = quatDot + rates[4:7, 0] = posDot + rates[7:13, 0] = PDot + return rates + +def fWrapper(y, t, parameters): + y = mat(y).T + dy = rateOfChange(y, t, parameters) + return dy.T.A1 + +def SCIPY(endtime, dt, state, parameters): + times = np.arange(0.0, endtime, dt) + y0 = state.T.A1 + res = scipy.integrate.odeint(fWrapper, y0, times, args=(parameters,)) + states = [] + res = res.T + r,c = res.shape + for ci in range(0,c): + states.append(mat(res[:,ci]).T) + return states + +def RK4(endtime, dt, state, parameters): + t = 0.0 + states = [] + while t < endtime: + newstate = mat (zeros( state.shape )) # Create a new object + + #### Runge Kutta integration: + k1 = rateOfChange(state, t, parameters) + k2 = rateOfChange(state + 0.5*dt*k1, t, parameters) + k3 = rateOfChange(state + 0.5*dt*k1, t, parameters) + k4 = rateOfChange(state + dt*k3, t, parameters) + newstate = state + (dt/6.0)*(k1+2*k2+2*k3+k4) + + # Normalize quat: + newstate[0:4, 0] = normalizeQuaternion(newstate[0:4, 0]) + states.append(newstate) + + state = newstate + + t += dt + print state[6,0], t, ' ', (t/endtime)*100.0, '%' + return states + + +def simulate(endtime, dt): + PInitial = mat( zeros((6,1)) ) + posInitial = mat(array([-1.2, -1.3, 2.8])).T + quatInitial = rotateAbout(mat(array([0.2, 1.0, 0.4])).T, 0.5) + # Parameters: + gravity = mat( array([0,0,-9.81]) ).T + massI = mat(eye(6)) * 1.01 # Mass matrix + parameters = (massI, gravity) + + # The state vector! + state = mat(zeros( (13,1) )) + state[0:4, 0] = quatInitial + state[4:7, 0] = posInitial + state[7:13, 0] = PInitial + + return SCIPY(endtime, dt, state, parameters) + +class W(QGLWidget): + time = 0.0 + index = 0 + def __init__(self, states, dt, parent=None): + super(W, self).__init__(parent) + self.firstRun = True + self.savePNGS = False + self.dt = dt + self.resize(500,500) + self.states = states + self.UP() + t = QTimer(self) + t.timeout.connect(self.UP) + t.start(self.dt*1000.0) + + def UP(self): + self.time += self.dt + if self.states: + state = self.states[self.index] + self.Dicequat = state[0:4, 0] + self.Dicepos = state[4:7, 0] + self.index += 1 + if self.index == len(self.states): + self.index = 0 + self.firstRun = False + # Paint: + self.update() + if self.firstRun: + # Create png images for the movie: + if self.savePNGS: + pm = self.renderPixmap() + pm.save('image'+str(self.index)+'.png') + + def initializeGL(self): + glClearColor(0.0, 0.5, 0.0, 1.0) + glEnable(GL_DEPTH_TEST) + glDepthFunc(GL_LESS) + glShadeModel(GL_SMOOTH) + def resizeGL(self,w,h): + glViewport(0, 0, w, h) + glMatrixMode(GL_PROJECTION) + glLoadIdentity() + gluPerspective(45.0, float(w)/float(h), 0.1, 100.0) + glMatrixMode(GL_MODELVIEW) + def paintGL(self): + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) # Clear buffers + glLoadIdentity() # Reset The View + + glLoadIdentity() + glTranslatef(0.0,-2.0,-10.0) # Move Left And Into The Screen + glRotatef(-90.0, 1.0, 0.0, 0.0) + drawFloor(2.0, 0.0) + drawAxis() + self.renderText(1.0, 0.0, 0.0, 'X') + self.renderText(0.0, 1.0, 0.0, 'Y') + self.renderText(0.0, 0.0, 1.0, 'Z') + + self.renderText(0.0,0.0,1.2,str(self.time)) + + x,y,z = self.Dicepos.A1 + R = quatToR(self.Dicequat) + + glTranslatef(x, y, z) + # Trick to rotate the openGL matrix: + r = R.A1 + rotR = (r[0], r[3], r[6], 0.0, + r[1], r[4], r[7], 0.0, + r[2], r[5], r[8], 0.0, + 0.0, 0.0, 0.0, 1.0) + glMultMatrixd(rotR) + + drawCube(0.6) + +et = 20.0 +dt = 0.04 +print 'starting integration... endtime =', et, ' stepsize =', dt +t0 = time.time() +states = simulate(et, dt) +t1 = time.time() +print 'That was heavy, it took me ', t1-t0, ' seconds!' +app = QApplication(sys.argv) +w = W(states, dt) +w.show() +sys.exit(app.exec_()) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/codeeditor.py Mon Dec 24 13:32:54 2012 +0100 @@ -0,0 +1,140 @@ +from PyQt4.QtCore import * +from PyQt4.QtGui import * +import compiler.lexer +import os.path + +class MySyntaxHighlighter(QSyntaxHighlighter): + def __init__(self, parent=None): + super(MySyntaxHighlighter, self).__init__(parent) + # Syntax highlighting: + self.rules = [] + fmt = QTextCharFormat() + fmt.setForeground(Qt.darkBlue) + fmt.setFontWeight(QFont.Bold) + for kw in compiler.lexer.keywords: + pattern = '\\b'+kw+'\\b' + self.rules.append( (pattern, fmt) ) + + # Comments: + fmt = QTextCharFormat() + fmt.setForeground(Qt.gray) + fmt.setFontItalic(True) + pattern = '\{.*\}' + self.rules.append( (pattern, fmt) ) + + # Procedure: + fmt = QTextCharFormat() + fmt.setForeground(Qt.blue) + fmt.setFontItalic(True) + #pattern = '(?<=procedure )[A-Za-z]' + # TODO lookbehind does not work, think something else + #self.rules.append( (pattern, fmt) ) + + def highlightBlock(self, text): + for pattern, fmt in self.rules: + expression = QRegExp(pattern) + index = expression.indexIn(text) + while index >= 0: + length = expression.matchedLength() + self.setFormat(index, length, fmt) + index = expression.indexIn(text, index + length) + +class LineNumberArea(QWidget): + def __init__(self, codeedit): + super(LineNumberArea, self).__init__(codeedit) + self.codeedit = codeedit + # TODO: display error in this: self.setToolTip('hello world') + def sizeHint(self): + return QSize(self.codeedit.lineNumberAreaWidth(), 0) + def paintEvent(self, ev): + self.codeedit.lineNumberAreaPaintEvent(ev) + +class CodeEdit(QPlainTextEdit): + def __init__(self, parent=None): + super(CodeEdit, self).__init__(parent) + # members: + self.isUntitled = True + self.filename = None + self.setFont(QFont('Courier')) + self.lineNumberArea = LineNumberArea(self) + + self.blockCountChanged.connect(self.updateLineNumberAreaWidth) + self.updateRequest.connect(self.updateLineNumberArea) + + # Syntax highlighter: + self.highlighter = MySyntaxHighlighter(self.document()) + + def setFileName(self, filename): + self.filename = filename + self.isUntitled = False + self.setWindowTitle(filename) + def setSource(self, source): + self.setPlainText(source) + + def save(self): + pass + def saveAs(self): + pass + + def saveFile(self): + if self.isUntitled: + self.saveAs() + else: + source = str(self.toPlainText()) + f = open(self.filename, 'w') + f.write(source) + f.close() + + def highlightErrorLocation(self, row, col): + tc = QTextCursor(self.document()) + tc.clearSelection() + tc.movePosition(tc.Down, tc.MoveAnchor, row - 1) + tc.movePosition(tc.Right, tc.MoveAnchor, col - 1) + tc.movePosition(tc.NextCharacter, tc.KeepAnchor) # Select 1 character + selection = QTextEdit.ExtraSelection() + lineColor = QColor(Qt.red).lighter(160) + selection.format.setBackground(lineColor) + #selection.format.setProperty(QTextFormat.FullWidthSelection, True) + selection.cursor = tc + self.setExtraSelections( [ selection ] ) + def clearErrors(self): + self.setExtraSelections( [ ] ) + + def lineNumberAreaWidth(self): + digits = 1 + mx = max(1, self.blockCount()) + while mx >= 10: + mx = mx / 10 + digits += 1 + space = 3 + self.fontMetrics().width('8') * digits + return space + def lineNumberAreaPaintEvent(self, ev): + painter = QPainter(self.lineNumberArea) + painter.fillRect(ev.rect(), Qt.lightGray) + block = self.firstVisibleBlock() + blockNumber = block.blockNumber() + top = self.blockBoundingGeometry(block).translated(self.contentOffset()).top() + bottom = top + self.blockBoundingRect(block).height() + while block.isValid() and top <= ev.rect().bottom(): + if block.isVisible() and bottom >= ev.rect().top(): + num = str(blockNumber + 1) + painter.setPen(Qt.black) + painter.drawText(0, top, self.lineNumberArea.width(), self.fontMetrics().height(), Qt.AlignRight, num) + block = block.next() + top = bottom + bottom = top + self.blockBoundingRect(block).height() + blockNumber += 1 + def resizeEvent(self, ev): + super(CodeEdit, self).resizeEvent(ev) + cr = self.contentsRect() + self.lineNumberArea.setGeometry(QRect(cr.left(), cr.top(), self.lineNumberAreaWidth(), cr.height() )) + def updateLineNumberAreaWidth(self, newBlockCount): + self.setViewportMargins(self.lineNumberAreaWidth(), 0, 0, 0) + def updateLineNumberArea(self, rect, dy): + if dy > 0: + self.lineNumberArea.scroll(0, dy) + else: + self.lineNumberArea.update(0, rect.y(), self.lineNumberArea.width(), rect.height()) + if rect.contains(self.viewport().rect()): + self.updateLineNumberAreaWidth(0) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/integration_methods.py Mon Dec 24 13:32:54 2012 +0100 @@ -0,0 +1,205 @@ +import sys +from random import random +import numpy as np +import pylab +import time +import scipy.integrate + +""" + This script compares different methods to solve an ODE. + + The solved system consists of: + y_dot = f(y, t, p) + y0 + + where + y: the state vector of the system + y0: the initial state of the system + y_dot: the derivatives of the state vector + f: the function that returns the derivatives. + t: the time instant + p: the parameters of the system +""" + +def func(y, t, parameters): + """ + The example system is a bouncing point mass on a hunt-cossley surface. + State vector: + y[0] = height + y[1] = vertical momentum + + y_dot[0] = velocity + y_dot[1] = force + """ + # Step 1: extract the states into variables: + h = y[0] + p = y[1] + mass = parameters[0] + Kp_floor = parameters[1] + Kd_floor = parameters[2] + + # Step 2: Calculate model: + v = p / mass + + # Ground reaction force: + if h > 0: + F = 0.0 + else: + F = -h * Kp_floor + + # Gravity: + F += -9.81 * mass + y_dot = np.array([v, F], dtype=np.float) + return y_dot + +parameters = np.array([3, 50000, 3000], dtype=np.float) +y0 = np.array([1, 0], dtype=np.float) # Start at 1 meter, with 0 momentum + +class SCIPY: + def __init__(self, dt): + self.dt = dt + def name(self): + return "scipy.integrate" + def integrate(self, f, y0, parameters, endtime): + timevector = np.arange(0.0, endtime, self.dt) + states = scipy.integrate.odeint(f, y0, timevector, args=(parameters,)) + return timevector, states + +class RK4: + def __init__(self, dt): + self.dt = dt + + def name(self): + return "RK4" + def integrate(self, f, y0, parameters, endtime): + print(y0, parameters) + t = 0.0 + states = [] + timevector = [] + state = y0 + states.append(state) + timevector.append(t) + while t < endtime: + #### Runge Kutta integration: + k1 = func(state, t, parameters) + k2 = func(state + 0.5*self.dt*k1, t, parameters) + k3 = func(state + 0.5*self.dt*k1, t, parameters) + k4 = func(state + self.dt*k3, t, parameters) + newstate = state + (self.dt/6.0)*(k1+2*k2+2*k3+k4) + + t += self.dt + state = newstate + + states.append(state) + timevector.append(t) + + return np.array(timevector, dtype=np.float), np.array(states, dtype=np.float) + + +class BDFfixed: + def __init__(self, order): + self.order = order + def name(self): + return "BDF" + def integrate(self, f, y0, parameters, endtime): + timevector = [] + states = [] + timevector.append(0.0) + states.append(y0) + t = 0.0 + h = 0.01 + state = y0 + while t < endtime: + if len(states) < 4: + k1 = func(state, t, parameters) + k2 = func(state + 0.5*h*k1, t, parameters) + k3 = func(state + 0.5*h*k1, t, parameters) + k4 = func(state + h*k3, t, parameters) + newstate = state + (h/6.0)*(k1+2*k2+2*k3+k4) + state = newstate + t += h + timevector.append(t) + states.append(state) + else: + t += h + if np.max(state) > 10: + break + hyn = (11.0/6.0)*states[-1] + (-3)*states[-2] + (3/2)*states[-3] + (-1/3)*states[-4] + state = state + hyn + timevector.append(t) + states.append(state) + + return np.array(timevector, dtype=np.float), np.array(states, dtype=np.float) + +class WODE: + def __init__(self): + pass + def integrate(self, f, y0, parameters, endtime): + timevector = [] + states = [] + timevector.append(0.0) + states.append(y0) + return np.array(timevector, dtype=np.float), np.array(states, dtype=np.float) + +class DOPRI45: + def __init__(self): + pass + def name(self): + return "DOPRI" + def integrate(self, f, y0, parameters, endtime): + timevector = [] + states = [] + timevector.append(0.0) + states.append(y0) + y = y0 + t = 0.0 + + h = 0.01 + iteration = 0 + while t < endtime: + k1 = h*f(y, t, parameters) + k2 = h*f(y + (1.0/5.0)*k1, t + (1.0/5.0)*h, parameters) + k3 = h*f(y + (3.0/40.0)*k1 + (9.0/40.0)*k2, t + (3.0/10.0)*h, parameters) + k4 = h*f(y + (44.0/45.0)*k1 - (56.0/15.0)*k2 + (32.0/9.0)*k3, t + (4.0/5.0)*h, parameters) + k5 = h*f(y + (19372.0/6561.0)*k1 - (25360.0/2187.0)*k2 + (64448.0/6561.0)*k3 - (212.0/729.0)*k4, t + (8.0/9.0)*h, parameters) + k6 = h*f(y + (9017.0/3168.0)*k1 - (355.0/33.0)*k2 - (46732.0/5247.0)*k3 + (49.0/176.0)*k4 - (5103.0/18656.0)*k5, t + h, parameters) + k7 = h*f(y + (35.0/384.0)*k1 + (500.0/1113.0)*k3 + (125.0/192.0)*k4 - (2187.0/6784.0)*k5 + (11.0/84.0)*k6, t + h, parameters) + + + y_next = y + (35.0/384.0)*k1 + (500.0/1113.0)*k3 + (125.0/192.0)*k4 - (2187.0/6784.0)*k5 + (11.0/84.0)*k6 + z_next = y + (5179.0/57600.0)*k1 + (7571.0/16695.0)*k3 + (393.0/640.0)*k4 - (92097.0/339200.0)*k5 + (187.0/2100.0)*k6 + (1.0/40.0)*k7 + error = np.linalg.norm(y_next - z_next) + #print(error, t) + eps = 0.01 + if error < eps: + y = y_next + t = t + h + timevector.append(t) + states.append(y) + + s = ((eps*h)/(2*error)) ** (1.0/5.0) + h = s * h + if h < 1e-5: + h = 1e-5 + iteration += 1 + if iteration > 10000: + break + + return np.array(timevector, dtype=np.float), np.array(states, dtype=np.float) + +rk4 = RK4(0.01) +sp = SCIPY(0.01) +wode = BDFfixed(4) +dopri = DOPRI45() +methods = [rk4, sp, dopri] + +pylab.figure() +pylab.hold(True) +for method in methods: + t, s = method.integrate(func, y0, parameters, 10) + print(t.shape, s.shape) + pylab.plot(t, s[:,0], 'x-', label=method.name()) + +pylab.legend() +pylab.show() +
--- a/python/libs/project.py Mon Dec 24 13:30:12 2012 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,34 +0,0 @@ -""" - source project - contains: - - modules - - primitives like functions, types and variables - - other modules -""" - -import json - -class Project: - def __init__(self, filename=None): - self.name = "" - self.modules = [] - self.elements = [] - self.settings = {} - if filename: - self.load(filename) - - def load(self, filename): - """ Load the project from file """ - self.filename = filename - - with open(self.filename, 'r') as f: - d = json.load(f) - self.elements = d['elements'] - - def save(self): - if self.filename: - d = {'elements': self.elements} - print(d) - with open(self.filename, 'w') as f: - json.dump(d, f) -
--- a/python/libs/widgets/__init__.py Mon Dec 24 13:30:12 2012 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -# File to make this directory a package. - -from .codeeditor import CodeEdit -from .astviewer import AstViewer -
--- a/python/libs/widgets/astviewer.py Mon Dec 24 13:30:12 2012 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,36 +0,0 @@ -from PyQt4.QtCore import * -from PyQt4.QtGui import * - -def astToNamedElement(astNode, parentNode): - """ Helper to convert and AST tree to NamedElement tree: """ - item = QStandardItem(str(astNode)) - item.setData(astNode) - parentNode.appendRow(item) - for c in astNode.getChildren(): - astToNamedElement(c, item) - -# The actual widget: -class AstViewer(QTreeView): - sigNodeSelected = pyqtSignal(object) - def __init__(self, parent=None): - super(AstViewer, self).__init__(parent) - self.setHeaderHidden(True) - self.clicked.connect(self.selectHandler) - - def setAst(self, ast): - """ Create a new model and add all ast elements to it """ - model = QStandardItemModel() - if ast: - astToNamedElement(ast, model.invisibleRootItem()) - self.setModel( model ) - self.expandAll() - - def selectHandler(self, index): - if not index.isValid(): - return - model = self.model() - item = model.itemFromIndex(index) - node = item.data() - self.sigNodeSelected.emit(node) - -
--- a/python/libs/widgets/codeeditor.py Mon Dec 24 13:30:12 2012 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,140 +0,0 @@ -from PyQt4.QtCore import * -from PyQt4.QtGui import * -import compiler.lexer -import os.path - -class MySyntaxHighlighter(QSyntaxHighlighter): - def __init__(self, parent=None): - super(MySyntaxHighlighter, self).__init__(parent) - # Syntax highlighting: - self.rules = [] - fmt = QTextCharFormat() - fmt.setForeground(Qt.darkBlue) - fmt.setFontWeight(QFont.Bold) - for kw in compiler.lexer.keywords: - pattern = '\\b'+kw+'\\b' - self.rules.append( (pattern, fmt) ) - - # Comments: - fmt = QTextCharFormat() - fmt.setForeground(Qt.gray) - fmt.setFontItalic(True) - pattern = '\{.*\}' - self.rules.append( (pattern, fmt) ) - - # Procedure: - fmt = QTextCharFormat() - fmt.setForeground(Qt.blue) - fmt.setFontItalic(True) - #pattern = '(?<=procedure )[A-Za-z]' - # TODO lookbehind does not work, think something else - #self.rules.append( (pattern, fmt) ) - - def highlightBlock(self, text): - for pattern, fmt in self.rules: - expression = QRegExp(pattern) - index = expression.indexIn(text) - while index >= 0: - length = expression.matchedLength() - self.setFormat(index, length, fmt) - index = expression.indexIn(text, index + length) - -class LineNumberArea(QWidget): - def __init__(self, codeedit): - super(LineNumberArea, self).__init__(codeedit) - self.codeedit = codeedit - # TODO: display error in this: self.setToolTip('hello world') - def sizeHint(self): - return QSize(self.codeedit.lineNumberAreaWidth(), 0) - def paintEvent(self, ev): - self.codeedit.lineNumberAreaPaintEvent(ev) - -class CodeEdit(QPlainTextEdit): - def __init__(self, parent=None): - super(CodeEdit, self).__init__(parent) - # members: - self.isUntitled = True - self.filename = None - self.setFont(QFont('Courier')) - self.lineNumberArea = LineNumberArea(self) - - self.blockCountChanged.connect(self.updateLineNumberAreaWidth) - self.updateRequest.connect(self.updateLineNumberArea) - - # Syntax highlighter: - self.highlighter = MySyntaxHighlighter(self.document()) - - def setFileName(self, filename): - self.filename = filename - self.isUntitled = False - self.setWindowTitle(filename) - def setSource(self, source): - self.setPlainText(source) - - def save(self): - pass - def saveAs(self): - pass - - def saveFile(self): - if self.isUntitled: - self.saveAs() - else: - source = str(self.toPlainText()) - f = open(self.filename, 'w') - f.write(source) - f.close() - - def highlightErrorLocation(self, row, col): - tc = QTextCursor(self.document()) - tc.clearSelection() - tc.movePosition(tc.Down, tc.MoveAnchor, row - 1) - tc.movePosition(tc.Right, tc.MoveAnchor, col - 1) - tc.movePosition(tc.NextCharacter, tc.KeepAnchor) # Select 1 character - selection = QTextEdit.ExtraSelection() - lineColor = QColor(Qt.red).lighter(160) - selection.format.setBackground(lineColor) - #selection.format.setProperty(QTextFormat.FullWidthSelection, True) - selection.cursor = tc - self.setExtraSelections( [ selection ] ) - def clearErrors(self): - self.setExtraSelections( [ ] ) - - def lineNumberAreaWidth(self): - digits = 1 - mx = max(1, self.blockCount()) - while mx >= 10: - mx = mx / 10 - digits += 1 - space = 3 + self.fontMetrics().width('8') * digits - return space - def lineNumberAreaPaintEvent(self, ev): - painter = QPainter(self.lineNumberArea) - painter.fillRect(ev.rect(), Qt.lightGray) - block = self.firstVisibleBlock() - blockNumber = block.blockNumber() - top = self.blockBoundingGeometry(block).translated(self.contentOffset()).top() - bottom = top + self.blockBoundingRect(block).height() - while block.isValid() and top <= ev.rect().bottom(): - if block.isVisible() and bottom >= ev.rect().top(): - num = str(blockNumber + 1) - painter.setPen(Qt.black) - painter.drawText(0, top, self.lineNumberArea.width(), self.fontMetrics().height(), Qt.AlignRight, num) - block = block.next() - top = bottom - bottom = top + self.blockBoundingRect(block).height() - blockNumber += 1 - def resizeEvent(self, ev): - super(CodeEdit, self).resizeEvent(ev) - cr = self.contentsRect() - self.lineNumberArea.setGeometry(QRect(cr.left(), cr.top(), self.lineNumberAreaWidth(), cr.height() )) - def updateLineNumberAreaWidth(self, newBlockCount): - self.setViewportMargins(self.lineNumberAreaWidth(), 0, 0, 0) - def updateLineNumberArea(self, rect, dy): - if dy > 0: - self.lineNumberArea.scroll(0, dy) - else: - self.lineNumberArea.update(0, rect.y(), self.lineNumberArea.width(), rect.height()) - if rect.contains(self.viewport().rect()): - self.updateLineNumberAreaWidth(0) -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/project.py Mon Dec 24 13:32:54 2012 +0100 @@ -0,0 +1,34 @@ +""" + source project + contains: + - modules + - primitives like functions, types and variables + - other modules +""" + +import json + +class Project: + def __init__(self, filename=None): + self.name = "" + self.modules = [] + self.elements = [] + self.settings = {} + if filename: + self.load(filename) + + def load(self, filename): + """ Load the project from file """ + self.filename = filename + + with open(self.filename, 'r') as f: + d = json.load(f) + self.elements = d['elements'] + + def save(self): + if self.filename: + d = {'elements': self.elements} + print(d) + with open(self.filename, 'w') as f: + json.dump(d, f) +