Installing CircuitPython on Arduino MKR WIFI

Part I: Installing CircuitPython

In this part, we will be flashing two bootloaders sequentially onto MKR.

1. Download and install the latest Arduino IDE https://www.arduino.cc/en/Main/Software
2. Open the IDE. Go to Tools -> Boards -> Boards Manager. In Boards Manager, search "Arduino SAMD Boards (32-bits AMR Cortex-M0+)" and then click Install. After installation, select Tools -> Boards -> MKR WIFI 1010
3. Plug-in MKR WIFI into the computer via a USB cable, and select Tools -> Port -> MKR WIFI or the only USB device when MKR is plugged in.
4. Download latest INO for MKR Zero's bootloader from https://github.com/adafruit/uf2-samdx1/releases to your computer, e.g., https://github.com/adafruit/uf2-samdx1/releases/download/v3.9.0/update-bootloader-mkrzero-v3.9.0.ino
5. Open the INO file just downloaded and upload it to MKR. You should see something like this on your Arduino IDE:

Sketch uses 21112 bytes (8%) of program storage space. Maximum is 262144 bytes.
Global variables use 3464 bytes of dynamic memory.
Atmel SMART device 0x10010005 found
Device : ATSAMD21G18A
Chip ID : 10010005
Version : v2.0 [Arduino:XYZ] Mar 19 2018 09:45:14
Address : 8192
Pages : 3968
Page Size : 64 bytes
Total Size : 248KB
Planes : 1
Lock Regions : 16
Locked : none
Security : false
Boot Flash : true
BOD : true
BOR : true
Arduino : FAST_CHIP_ERASE
Arduino : FAST_MULTI_PAGE_WRITE
Arduino : CAN_CHECKSUM_MEMORY_BUFFER
Erase flash
done in 0.823 seconds

Write 21112 bytes to flash (330 pages)

[===== ] 19% (64/330 pages)
[=========== ] 38% (128/330 pages)
[================= ] 58% (192/330 pages)
[======================= ] 77% (256/330 pages)
[============================= ] 96% (320/330 pages)
[==============================] 100% (330/330 pages)
done in 0.128 seconds

Verify 21112 bytes of flash with checksum.
Verify successful
done in 0.018 seconds
CPU reset.


6. MKR will restart. If not, press the reset button. A USB storage device called MKRZEROBOOT shall pop-up on your computer now. If a file called

CURRENT.UF2

is there, then you are good.

7. Download the CircuitPython image for MKR from https://circuitpython.org/board/arduino_mkrzero/, the UF2 file. Then copy it to the MKRZEROBOOT device. Then the MKR will restart and the new USB storage device's name will become CIRCUITPY. If it doesn't restart automatically, press the reset button.


Part II: Play around with CircuitPython

1. Download and install Mu editor. For Windows and Mac, go to https://codewith.mu/en/download . For Linux, or Python on other platforms, you can simply use pip: pip3 install mu-editor --user Then start the Mu editor.

2. Upon initial start, select "Adafruit CircuitPython". Then you see an blank editor. Now, click the "Serial" button at the top of the window. And you shall see the "Adafruit CircuitPython REPL" box. Click in it, and press any key to activate the Python console. Then, you can enjoy the Python shell.

3. To run a script, in the Mu editor, paste the code below and save it as a file on your CIRCUITPY storage device and make sure it is the only file there. After saving, the MKR will restart on its own -- if not, press the reset button. And you shall see the orange LED blinking.

import board
import digitalio
import time

led = digitalio.DigitalInOut(board.D6)
led.direction = digitalio.Direction.OUTPUT

while True:
print("Hello, CircuitPython!")
led.value = True
time.sleep(1)
led.value = False
time.sleep(1)


Note that the interactive mode Python shell and the scripting mode (i.e., running code from a script file) cannot coexist.

Jupyter server OSError: [Errno 99] Cannot assign requested address

If you run into this error:

OSError: [Errno 99] Cannot assign requested address

a likely fix is going to /etc/hosts and pointing localhost only to ipv4 loopback (127.0.0.1), not to ipv6's (::1).

For example, on my computer, /etc/hosts looks like this:

127.0.0.1       localhost

# The following lines are desirable for IPv6 capable hosts
::1     localhost       ip6-localhost   ip6-loopback
ff02::1 ip6-allnodes
ff02::2 ip6-allrouters

xxx.xxx.xxx.xxx  tinyai  tinyai

Just change it to

127.0.0.1       localhost

# The following lines are desirable for IPv6 capable hosts
::1           ip6-localhost   ip6-loopback
ff02::1 ip6-allnodes
ff02::2 ip6-allrouters

xxx.xxx.xxx.xxx  tinyai  tinyai

Note the change at the ::1 line.

For more details, read https://github.com/ipython/ipython/issues/6193#issuecomment-466082026

Python lists as local variables in recursion

Today I ran into a really interesting problem. Let's use Leetcode's Problem 77 combination as an example. I wanted write a function that prints all possible combinations (no order, unlike permutation).

So here is my function definition and call:

def helper(n, k, start, L):
    # L = list (L)  # make a copy 
    if len(L) == k:# all k elements have been drawn
        print L
    else:
        for i in range(start, n+1): 
            L.append(i)
            print "L before recursion", l
            helper(n,k,i+1, L) # down to next level
            print "L after recursion", l
            L = L[:-1] 

helper(3,2,1, []) # picking 2 from 3

But this was what I got:

L before recursion [1]
L before recursion [1, 2]
[1, 2]
L after recursion [1, 2]
L before recursion [1, 3]
[1, 3]
L after recursion [1, 3]   
L after recursion [1, 2]  --- l is supposed to be [1] here as quitting the inner recursion 
L before recursion [1, 2]
[1, 2]
L after recursion [1, 2]
L before recursion [1, 3]
[1, 3]
L after recursion [1, 3]

As you can see, L never came back to [1]. It became (at the 8th line of the print-out) [1,2] after executing helper(3,2,2, [1]). Further, there are unmatching numbers of "before" and "after" which made no sense to me.

It took me a couple of hours to think until I realized that maybe the reference of a list was the culprit. So I uncommented the "make a copy" line, and then things worked as expected:

L before recursion [1]
L before recursion [1, 2]
[1, 2]
L after recursion [1, 2]
L before recursion [1, 3]
[1, 3]
L after recursion [1, 3]
L after recursion [1]
L before recursion [2]
L before recursion [2, 3]
[2, 3]
L after recursion [2, 3]
L after recursion [2]
L before recursion [3]
L after recursion [3]

But I still had no clue why a copy fixed it, and even the in-and-out of recursions.
If you know why, let me know.

PS: An alternative fix is to let the function to return L[:-1]. In this way, the old l can be restored after a recursive call.

Named Entity Recognition using SpaCy in 5 minutes

Recently, I am looking it SpaCy, a startup and an NLP toolkit. It is fabulous on its speed. Today, I just gave it a try on NER. Just a few lines (as in iPython):

In [1]: import spacy.en
In [2]: parser = spacy.en.English()
In [12]: ParsedSentence = parser(u"alphabet is a new startup specializing in eating their own words on leaving china to fight for information freedom")

In [13]: print ParsedSentence.ents
()

In [14]: ParsedSentence = parser(u"Alphabet is a new startup specializing in eating their own words on leaving China to fight for information freedom")

In [15]: print ParsedSentence.ents
(Alphabet, China)

In [16]: for Entity in  ParsedSentence.ents:    
   ....:     print Entity.label, Entity.label_, ' '.join(t.orth_ for t in Entity)
   ....:     
349 ORG Alphabet
350 GPE China

I used only default settings. Apparently, the NER of SpaCy is very sensitive to case of words.

MinMaxScaler() of sklearn.preprocessing

This is just a note demonstrating how to scale and de-scale data using MinMaxScaler() of sklearn.preprocessing.

In [4]: import sklearn.preprocessing

In [5]: X=sklearn.preprocessing.MinMaxScaler()

In [7]: X.fit_transform([1.,2.,3.,4.])
Out[7]: array([ 0.        ,  0.33333333,  0.66666667,  1.        ])

In [10]: X.scale_
Out[10]: 0.33333333333333331

In [14]: X.min_
Out[14]: -0.33333333333333331

In [17]: 4.*X.scale_+X.min_  # from non-scaled value to scaled value, a.k.a., transform
Out[17]: 1.0

In [19]: (1-X.min_)/X.scale_ # from scaled value back to scaled_value, a.k.a., inverse_transform
Out[19]: 4.0

References:

[1] http://stackoverflow.com/questions/24724717/scikit-learn-minmaxscaler-produces-slightly-different-results-than-a-numpy-imple

Guide to CVXOPT's quadprog() for row-major and/or MATLAB-speaking minds

I am using CVXOPT now. My mind is very row-major, because I mostly program in C and Python. And my mind is also MATLAB-style, because my friend send me a code in MATLAB. Now I need to convert his code to CVXOPT in Python. Here are some thing you wanna pay attention.

1. Matrixes in CVXOPT are column-major

By the code:

G=cvxopt.matrix([[1.,1],[-1,2],[2,1],[-1,0],[0,-1]])

I thought I created a 5-by-2 matrix (5 horizontal rows and 2 vertical columns). However, I got a 2-by-5 matrix.

In [72]: G
Out[72]: <2x5 matrix, tc='d'>

In [73]: print G
[ 1.00e+00 -1.00e+00  2.00e+00 -1.00e+00  0.00e+00]
[ 1.00e+00  2.00e+00  1.00e+00  0.00e+00 -1.00e+00]

2. The matrixes must be of type float.

Make sure your matrixes are of type float, even though your coefficients are integers.

Earlier I created a matrix:

In [53]: q=cvxopt.matrix([[-2],[-6]])

Later I got the error like

TypeError: 'q' must be a 'd' matrix with one column

This is because the matrix is not of type float:

In [54]: q
Out[54]: <1x2 matrix, tc='i'>

To fix, simply make at least one element float, e.g., "1." Then the type code will become d.

In [89]: q
Out[89]: <1x2 matrix, tc='d'>

3. Generating the matrixes needed for your optimization problem.

First of all, the cvxopt formulates quadratic programming so much simpler than MATLAB. Here is how cvxopt formulates the problem:

$$ \begin{array}{ll}
\mbox{minimize} & (1/2) x^TPx + q^T x \\
\mbox{subject to} & G x \preceq h, \\
& Ax = b.
\end{array}
$$

The boundaries for all variables can be defined as part of the matrix G - just use eyes. E.g., the lower boundary constraints

$$x_1 > 1, x_2>2, x_3> 3
$$
is just
$$
\begin{bmatrix}
-1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 0 & -1
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2 \\
x_ 3
\end{bmatrix}
< \begin{bmatrix} -1\\ -2\\ -3 \end{bmatrix} $$ Note that in either MATLAB's or cvxopt's formulation, the linear terms are always on the side of "smaller than or equal to". Hence, the minus one above. MATLAB formulations have two additional column arrays for upper and lower boundaries: $$ \begin{array}{ll} \mbox{minimize} & (1/2) x^THx + f^T x \\ \mbox{subject to} & A x \preceq b, \\ & A_{eq}x = b_{eq}, \\ & lb \preceq x \preceq ub. \end{array} $$ Apparently, there are lazy people who wants to have a straightforward conversion from MATLAB's formulation to cvxopt's. However, I am not sure whether that translation is out-dated because cvxopt changed their APIs. Here is an easier one.

Let the MATLAB version be

x = quadprog(H,f,A,b,Aeq,beq,lb,ub, ...)

This is how you do it in cvxopt:

import numpy
import cvxopt
import cvxopt.solvers
n = H.shape[1]   # n is the number of variables 

P = H
q = f
G = numpy.vstack([A, -numpy.eye(n), numpy.eye(n)])
h = numpy.vstack([b, -lb, ub])
A = Aeq
b = beq

sol = cvxopt.solvers.qp(cvxopt.matrix(P), cvxopt.matrix(q), cvxopt.matrix(G), cvxopt.matrix(h), cvxopt.matrix(A), cvxopt.matrix(b))
x = sol['x']

If you don't have lb nor/and ub, no need to have it/them in G and h.

Put things together

So, let's see an example. The example in MATLAB Optimization Toolbox can be solved in cvxopt as follows (in iPython shell):

In [81]: G=cvxopt.matrix([[1.,1],[-1,2],[2,1],[-1,0],[0,-1]])

In [82]: P=cvxopt.matrix([[1,-1],[-1,2.]])

In [83]: q=cvxopt.matrix([[-2,-6.]])

In [84]: h=cvxopt.matrix([2,2,3,0,0.])

In [85]: Solv=cvxopt.solvers.qp(P.T, q, G.T, h)
     pcost       dcost       gap    pres   dres
 0: -1.0389e+01 -8.2778e+00  2e+01  9e-01  1e+00
 1: -7.2856e+00 -9.9286e+00  3e+00  1e-16  4e-16
 2: -8.1161e+00 -8.6188e+00  5e-01  8e-17  2e-16
 3: -8.2068e+00 -8.2359e+00  3e-02  6e-17  3e-15
 4: -8.2220e+00 -8.2224e+00  3e-04  5e-17  1e-15
 5: -8.2222e+00 -8.2222e+00  3e-06  8e-17  6e-16
Optimal solution found.

In [86]: print Solv['x']
[ 6.67e-01]
[ 1.33e+00]

Reading Physiology Data in EDF+ Format Using EEGTools

For years, I always converted EDF+ format data into text files for processing. Recently, I was tired of this and decided to try directly read EDF+ data. I tried out a Python module EEGTools by Boris Reuderink and it turned out to be great. It's super easy to use. Just use these two lines:

import eegtools
Data = eegtools.io.load_edf('Sample.edf')

The structure of Data is like this

EDF(X=array([[ -14.,    9.,  -70., ...,  -50.,  -99.,  -93.],
       [ -35.,  -21.,  -87., ...,  -45.,  -98.,  -81.],
       [ -72.,  -62., -116., ...,  -32.,  -79.,  -52.],
       ..., 
       [ -37.,  -89.,  -63., ...,  -41.,  -16.,  -45.],
       [ -69., -120.,  -99., ...,  -46.,  -25.,  -55.],
       [ -44.,  -87.,  -57., ...,  -29.,   -4.,  -41.]], dtype=float32), sample_rate=160.0, chan_lab=['Fc5.', 'Fc3.', 'Fc1.', 'Fcz.', 'Fc2.', 'Fc4.', 'Fc6.', 'C5..', 'C3..', 'C1..', 'Cz..', 'C2..', 'C4..', 'C6..', 'Cp5.', 'Cp3.', 'Cp1.', 'Cpz.', 'Cp2.', 'Cp4.', 'Cp6.', 'Fp1.', 'Fpz.', 'Fp2.', 'Af7.', 'Af3.', 'Afz.', 'Af4.', 'Af8.', 'F7..', 'F5..', 'F3..', 'F1..', 'Fz..', 'F2..', 'F4..', 'F6..', 'F8..', 'Ft7.', 'Ft8.', 'T7..', 'T8..', 'T9..', 'T10.', 'Tp7.', 'Tp8.', 'P7..', 'P5..', 'P3..', 'P1..', 'Pz..', 'P2..', 'P4..', 'P6..', 'P8..', 'Po7.', 'Po3.', 'Poz.', 'Po4.', 'Po8.', 'O1..', 'Oz..', 'O2..', 'Iz..'], time=array([  0.00000000e+00,   6.25000000e-03,   1.25000000e-02, ...,
         1.22981250e+02,   1.22987500e+02,   1.22993750e+02]), annotations=[(0.0, 4.1, [u'T0']), (4.1, 4.1, [u'T1']), (8.2, 4.1, [u'T0']), (12.3, 4.1, [u'T2']), (16.4, 4.1, [u'T0']), (20.5, 4.1, [u'T2']), (24.6, 4.1, [u'T0']), (28.7, 4.1, [u'T1']), (32.8, 4.1, [u'T0']), (36.9, 4.1, [u'T2']), (41.0, 4.1, [u'T0']), (45.1, 4.1, [u'T1']), (49.2, 4.1, [u'T0']), (53.3, 4.1, [u'T2']), (57.4, 4.1, [u'T0']), (61.5, 4.1, [u'T1']), (65.6, 4.1, [u'T0']), (69.7, 4.1, [u'T1']), (73.8, 4.1, [u'T0']), (77.9, 4.1, [u'T2']), (82.0, 4.1, [u'T0']), (86.1, 4.1, [u'T1']), (90.2, 4.1, [u'T0']), (94.3, 4.1, [u'T2']), (98.4, 4.1, [u'T0']), (102.5, 4.1, [u'T1']), (106.6, 4.1, [u'T0']), (110.7, 4.1, [u'T2']), (114.8, 4.1, [u'T0']), (118.9, 4.1, [u'T1'])])

As you can see, the recordings are in the member X. The channels that member arrays of X correspond to are sequentially given in the member chan_lab. Sampling rate and time instants are given in sample_rate and time. Lastly, you can access the annotation on EEG from annotations.

To download or install, see EEGTools' GitHub page: http://borisreuderink.nl/ Read the instruction at the bottom and use PIP to install - easier than cloning and setup.py.

Signal Filtering using inverse FFT in Python

A straight forward way of doing signal filtering is zeroing out terms in inverse FFT result. I saw a good post online. However, it does not encapsulate into a function nor allow users to specify passing bands in terms of physical frequency. So I modified it into a function below.

def bandpass_ifft(X, Low_cutoff, High_cutoff, F_sample, M=None):
    """Bandpass filtering on a real signal using inverse FFT
    
    Inputs
    =======
    
    X: 1-D numpy array of floats, the real time domain signal (time series) to be filtered
    Low_cutoff: float, frequency components below this frequency will not pass the filter (physical frequency in unit of Hz)
    High_cutoff: float, frequency components above this frequency will not pass the filter (physical frequency in unit of Hz)
    F_sample: float, the sampling frequency of the signal (physical frequency in unit of Hz)    
    
    Notes
    =====
    1. The input signal must be real, not imaginary nor complex
    2. The Filtered_signal will have only half of original amplitude. Use abs() to restore. 
    3. In Numpy/Scipy, the frequencies goes from 0 to F_sample/2 and then from negative F_sample to 0. 
    
    """        
    
    import scipy, numpy
    if M == None: # if the number of points for FFT is not specified
        M = X.size # let M be the length of the time series
    Spectrum = scipy.fft(X, n=M) 
    [Low_cutoff, High_cutoff, F_sample] = map(float, [Low_cutoff, High_cutoff, F_sample])
    
    #Convert cutoff frequencies into points on spectrum
    [Low_point, High_point] = map(lambda F: F/F_sample * M /2, [Low_cutoff, High_cutoff])# the division by 2 is because the spectrum is symmetric 

    Filtered_spectrum = [Spectrum[i] if i >= Low_point and i <= High_point else 0.0 for i in xrange(M)] # Filtering
    Filtered_signal = scipy.ifft(Filtered_spectrum, n=M)  # Construct filtered signal 
    return Spectrum, Filtered_spectrum, Filtered_signal, Low_point, High_point

Here is an example showing its usage.

First of all, let's generate a signal of different frequencies.

import numpy
N = 400 # signal length of number
x = numpy.arange(0, N, 1) # generate the time ticks
Sines = [numpy.sin(x*n)*(1-n) for n in [.9, .75, .5, .25, .12, .03, 0.025]]# different frequency components
y = numpy.sum(Sines, axis=0)   # add them by column, low frequencies have higher amplitudes

import matplotlib.pyplot as plt
plt.plot(x, y) # visualize the data

The raw signal looks like this:



Then, let's assume the sampling frequency is 500Hz and we only want signals between 5 to 30 Hz to go thru the filter.

Low_cutoff, High_cutoff, F_sample = 5, 30, 500
Spectrum, Filtered_spectrum, Filtered_signal, Low_point, High_point = bandpass_ifft(y, Low_cutoff, High_cutoff, F_sample)

# Below is visualization
fig1 = plt.figure()
plt.stem(numpy.fft.fftfreq(N)*F_sample, abs(Spectrum))
plt.axis([0, F_sample/2, None, None]) # since the signal is real, the spectrum is symmetric|
plt.axvspan(Low_cutoff, High_cutoff,fc='g',alpha='.5') # The frequencies that we wanna get rid of

fig2 = plt.figure()
plt.plot(x, Filtered_signal)

In the first figure plotted, the green box shows frequencies that will be kept by the filter.



Here is the filtered signal. As you can see, high frequencies are gone. The signal looks very smooth.


Of course, this approach is very costly: FFT and iFFT are both $O(n\lg n)$ complexity algorithms. It is much more costly than convolving signal with filter terms (what people normally do in digital signal processing). Oh, and it cannot be done in real-time - you must wait until all data is sampled.

There is also a minor implementation issue here. Both Scipy and Numpy provides FFT/iFFT functions. But their speeds vary [1]. If you wanna go faster, try Python binding to some C/C++-based FFT libraries.

References:
1. http://stackoverflow.com/questions/6363154/what-is-the-difference-between-numpy-fft-and-scipy-fftpack

Stepper motor control with EasyDriver and Python on Linux desktop

My Mechanical Engineering colleagues recently asked me to build them a software-controlled stepper motor system - they wanna apply high-power pulse laser on the surface of rotating objects. Of course, the software must have GUI.

Hardware

I decided to control the stepper motor using a Linux PC with parallel ports (yes, people are still building motherboards with parallel ports), because developing software for desktop Linux is much easier. Oh, for tablets and smart phone fans, how do your devices talk to physical world? 

In searching for a proper driver chip/board for the stepper motor, I came across the handy EasyDriver board designed by Brian Schmalzhaus. I just need to provide one rising edge to the STEP pin of it, the stepper motor will rotate one step. I can control the direction of rotation and step size by configuring certain pins. It saves me the hassle of computing waveforms for controlling stepper motors. Check the references below to see how EasyDriver can simply your life.

Lastly, very important, let the EasyDriver board and the PC share the common ground. Otherwise the control signal from parallel port may not work (e.g., a high from parallel port is not high to the ground level of EasyDriver), or the parallel port might get burned.

Oh, one more thing, the EasyDriver board only drives bipolar stepper motors.

My circuit is prototyped on a breadboard like below:

The circuit takes power from the power supply of the desktop computer on which the control software will run. I used Pin 2 and Pin 4 of the parallel port to control the STEP and DIR on Easy Driver. MS1 and MS2 are grounded because I only use full step mode. I used two LEDs to show the signals from parallel ports.

Here is circuit in reality. As you can see, the circuit takes power from the old floppy disk power jack.

I used Saleae's Logic as a USB logic analyzer to check the waveforms. As you can see, for every rising edge on channel 4, the state of the first 4 channels change.

Software

At first I planned to use parapin, a C-based library for parallel port I/O. I used it when I was in college. After googling, I found out a Python module called PyParallel and it greatly simplifies software development.

I wrote a small piece of Python code below to control it.

"""This software is free software licensed under GNU GPL 3.0
Copyleft 2014 Forrest Sheng Bao http://fsbao.net

Dependencies: 
1. pyparallel: Can be installed on Ubuntu via
            sudo apt-get install python-parallel

Notes for Linux: 
1. make sure you execute
      sudo rmmod lp
           and
      sudo modprobe ppdev
  on your Unix shell first. 
2. Add yourself into lp user group on your system if you don't wanna use sudo all the time. 

Help the people under Communism tyranny in China, North Korea, Viet Nam and Cuba. 

"""

def rotate(Step, Dir, Step_Pin=5, Dir_Pin=2, Delay=0.25):
    """Rotate the motor by _Step_ steps in _Dir_ direction

    Input
    =======
    Step: integer
    Dir: string, 0 for clockwise and 1 for counterclockwise
    Step_Pin: int, the pin used to control stepping signal
    Dir_Pin: int, the pin used to control rotation direction
    Delay: float, the delay between each step. Do not set it 
           too small to exceed your motor's speed limit. 

    Notes
    =======
    The data pins in parallel ports are from 2 to 8. 

    """

    import parallel
    import time
    p = parallel.Parallel()

    p.setData(0)

    for i in range(int(Step)):
        p.setData(2 ** (Dir_Pin - 2) * Dir) # set or keep direction
        time.sleep(Delay)
        p.setData(2 ** (Dir_Pin - 2) * Dir + 2 ** (Step_Pin - 2)) # keep direction and rotate
        time.sleep(Delay)
    
#    p.setData(0)
    
if __name__ == "__main__" :
    import sys
    Step = int(sys.argv[1])
    Dir  = int(sys.argv[2])
    Delay = float(sys.argv[3])
    rotate(Step, Dir, Delay=Delay)

Note

1. The chip on EasyDriver can get really hot - up to 136 Celsius! Do not use it with a motor that requires lots of current and/or rotate the motor too frequently. Mine got burned out pretty soon. So I am ordering a motor with less torque but smaller current.
2. Do not send pulse to STEP pin of EasyDriver faster than the maximum speed of your motor. It will cause damage.

Questions and comments are all welcome!

References:

1. Easy Driver Example by Brain Schmalzhaus himself, http://www.schmalzhaus.com/EasyDriver/Examples/EasyDriverExamples.html
2. Easy Driver Example by Sparkfun, https://www.sparkfun.com/tutorials/400

An incompatible problem between python-vtk and pyvtk

I am using python-vtk (official Python binding of VTK) and pyvtk (only for VTK-format file I/O) these days. It's not a pleasant experience. It took me a few hours to debug my code and found out the problem at somewhere I never expected - the incompatibility of the libraries I am using. On top of that, one of them has compatibility issue with Python2.7.

I used python-vtk to read from VTK files (because I couldn't google out examples using pyvtk to do so). And then I used pyvtk to write into VTK files (because I couldn't google out examples using python-vtk to create scalar POINTDATA).

The problem is that in python-vtk, vtk.vtkDataSetReader().GetOutput().GetPoint() returns a 3-tuple, which is the X-, Y- and Z-coordinates of a POINT in DATASET POLYDATA block, whereas in pyvtk, pyvtk.PolyData() takes in coordinates as a list, not a 3-tuple.

I don't understand how this could cause problems, because lists and tuples are very similar in Python. After mandatory type conversion, the problem was gone.

pyvtk also has problem with Python2.7. I haven't debugged out the cause. But I just know it does not work with Python2.7.

PS: There is a great lack of resources on these two libraries. I use one library for reading VTK files and the other for writing VTK files, though each of them has the ability to do both. This "complimentary" combination is because I couldn't google out examples on using the opposite library for the opposite function. I really don't understand this. Without document, a software package is nothing.

I didn't find any documentation on pyvtk's website. For python-vtk, released with VTK (C++ mainly, along with Tcl, Python and Java bindings), I still didn't find docs for Python. I had to use iPython to learn.

Does anyone know why?

Reading VTK files in Python via python-vtk

by Forrest Sheng Bao http://fsbao.net

Update 2012-06-29: I just figured out a way to access polygons/cell using VTK's python binding. Check here: http://forrestbao.blogspot.com/2012/06/vtk-polygons-and-other-cells-as.html

Surprisingly, I noticed that there isn't a good document covering frequently-used functions of the Python module vtk (provided via python-vtk on Debian/Ubuntu Linux systems). So I decide to write a very simple one here, covering all functions that I have used.

Note:

• This tutorial is for Python 2.X though slight changes can make it work with Python 3.X.
• I assume you are familiar with VTK data format, thus you know what header, DATASET and POINT_DATA are.
• Notations like In [123] or Out[123] are prompts in iPython (an interactive Python shell). They show a line of code and its output/effect, respectively. They should NOT appear in your code. And, when you use other Python interpreter/shell, you may not see it.
• I did not use print function in iPython when showing something. But when you write a Python program, you need print to display.

Step 1: Set up the reader.

import vtk
reader = vtk.vtkDataSetReader()
reader.SetFileName("lh.sulc.fundi.from.pits.pial.vtk")
reader.ReadAllScalarsOn()  # Activate the reading of all scalars
reader.Update()

The reader is the top level object we use to access a VTK file. For example, you can know the header of a VTK file by

In [119]: reader.GetHeader()
Out[119]: 'Created by Mindboggle'

Step 2: Get data in DATASET block.

Again, reader is the top level. To get data in DATASET block, we need to call a method on reader.

data=reader.GetOutput()

We can know how many (not ``much'' here) data of each type are in the DATASET by the function GetNumberOf{VTK_Data_Type}(). For example, to know how may Points are there, we use:

In [11]: data.GetNumberOfPoints()
Out[11]: 128895

I used tab-completion function of iPython to find out all GetNumberOf{VTK_Data_Type}() for data:

data.GetNumberOfCells   data.GetNumberOfPolys
data.GetNumberOfLines   data.GetNumberOfStrips
data.GetNumberOfPieces  data.GetNumberOfVerts
data.GetNumberOfPoints  

After knowing the size of data of each type, we can access them. For example, to get the first point's coordinate, we can do this:

In [136]: data.GetPoint(0)
Out[136]: (-11.605026245117188, -97.47259521484375, 3.8222298622131348)

To access vertexes in the VERTICES block, it is a little different.

In [15]: data.GetNumberOfVerts()
Out[15]: 1L
In [18]: vt = data.GetVerts()
In [20]: vt.GetSize()
Out[20]: 56094L
In [23]: vt.GetData().GetValue(2)
Out[23]: 1091L
In [24]: vt.GetData().GetValue(3)
Out[24]: 1108L

Please note that some people, like myself, prefer to put all vertexes in one line, so the result of data.GetNumberOfVerts() and data.GetVerts().GetSize() look different. But it shouldn't bother you from using data.GetVerts().GetData().GetValue(x) to access them.

Step 3: get data in POINT_DATA/CELL_DATA block.

Data in POINT_DATA and CELL_DATA can be loaded via function GetPointData() and GetCellData(). For example,

d=data.GetPointData()

We will use loading a scalar array from POINT_DATA as an example below.
Vectors and tensors can be accessed in similar way by different functions.

All scalar arrays in a file can be viewed by:

In [140]: reader.GetNumberOfScalarsInFile() # get number of scalars
Out[140]: 5
In [141]: reader.GetScalarsNameInFile(1) # get scalar name string
Out[140]: curv

A scalar can be accessed by its name:

In [38]: array=d.GetArray('curv') # 'curv' is the scalar name
In [50]: array.GetValue(260957-260761+1)
Out[50]: 0.043181419372558594

In the end, you should see the following variables:

Variable   Type         Data/Info
---------------------------------
curv       vtkobject    vtkFloatArray (0x9e5cd38)<...>)\n  Array: 0x9f1a9a0\n\n
d          vtkobject    vtkPointData (0x9789b50)\<...>  PedigreeIds: (none)\n\n
data       vtkobject    vtkPolyData (0x9e6c650)\n<...>: 0\n  Ghost Level: 0\n\n
reader     vtkobject    vtkDataSetReader (0x9c431<...> InputStringLength: 0\n\n
vtk        module       < module 'vtk' from '/usr/<...>hon2.6/vtk/__init__.pyc'>

The end. Comments and questions are welcomed.

References:

• http://www.vtk.org/pipermail/vtkusers/2010-June/109492.html
• http://stackoverflow.com/questions/6684306/how-can-i-read-a-vtk-file-into-a-python-datastructure

Anything like MATLAB profiler for Python?

I just saw an interesting feature of MATLAB today. It's called Profiler in Parallel Computing Toolbox (http://www.mathworks.com/help/toolbox/distcomp/bra51qt-1.html). Basically, it's very similar to profile command (http://www.mathworks.com/help/techdoc/ref/profile.html).

You can know the time spending on each line of code. This is a very useful tool for scientific computing. I would like to know whether we have the similar stuff for Python.

To say "Hello world" in Python CGI Web Programming in 5 minutes

by Forrest Sheng Bao http://fsbao.net

I have lotta Python programs for bioinformatics research. I wanted to put them onto the web. I only developed Web apps in PHP before. And it seemed to be a big pain for porting a Python program to the web. But, I figured out in 5 minutes.

First, you need an Apache server on your Linux/Mac OS/Windows box. If you are gonna use data base, you need a database server. There are tons of blogs addressing these issues. So I won't be gossipy here. Suppose everything we mention below happens on your server-even web browsing.

Second, configure cgi in Apache. There are many ways to run a Python program on a web/http interface. I think CGI is the easiest. Assume you are on your own server and using all default settings. On Ubuntu Linux 9.04, the default configuration file for your default website is /etc/apache2/sites-available/default Open it, find the part for cgi directory, and make it like this

ScriptAlias /cgi-bin/ /var/www/cgi-bin/
<Directory "/var/www/cgi-bin">

AllowOverride None
Options +ExecCGI -MultiViews +SymLinksIfOwnerMatch
#Order allow,deny
Require all granted
Allow from all
AddHandler cgi-script .py              # tell Apache to handle every file with .py suffix as a cgi program
AddHandler default-handler .html .htm  # tell Apache to handle HTML files in regular way

</Directory>

The line

ScriptAlias /cgi-bin/ /var/www/cgi-bin/

specifies the path of cgi codes and the actual directory of your program. So when you type http://127.0.0.1/cgi-bin, the Apache will look into the directory /var/www/cgi-bin/ of your localhost.

Also, make sure that cgimodule of Apache is enabled. You can enable it by

sudo a2enmod cgi

Now restart your apache. On Ubuntu Linux, by default installation and configuration, it is

sudo /etc/init.d/apache2 restart

Third, write up a Hello, World! program.

#!/usr/bin/env python
print "Content-Type: text/html"
print
print """\
<html>
<head><title>First Python HTTP Programming </title></head>

<body>
<h2>Hello World!</h2>
</body>
</html>
"""

Now open http://127.0.0.1/cgi-bin/hello.py in your web browser and you shall see a hello world in it.

Reference: http://webpython.codepoint.net/cgi_tutorial

enumerate() array elements rather than indexing them in Python

by Forrest Sheng Bao http://fsbao.net

Recently, I am debugging two Python programs about graphs. I get very confused in many times, whether I am visiting an array element or the index of the array element. For example, I need to get the vertex number by indexing an array, and use that vertex number to index another array.

So, my friend Lay Yuan in Arizona State University told me a trick, using enumerate() in Python. This makes things much much more easier. I guess you are smart enough so this example can explain every thing to you

>>> for i,arr in enumerate([6,3,2,3,4]):
...     print i, arr
...
0 6
1 3
2 2
3 3
4 4

Below is an implementation of the Prim's algorithm. In this program, one step is to update the distance of current spanning tree to all neighbors of the newly added vertex. The neighbors of all vertexes are stored in an array neighbors. Without enumerate, I need these two lines to pick neighbors one by one:

for i in xrange(0,len(neighbors)):
w= neighbors[v][i];

With enumerate(), I just need one line to get the w:

for i, w in enumerate(neighbors[v]):

I don't need to struggle on the indexing problem any more.

# Python implementation of Prim's algorithm
# by Forrest Sheng Bao, Dept. of Computer Science, Texas Tech
# http://fsbao.net  forrest dot bao at gmail dot com
# This code implements the example on P. 571 of the famous
# "Introduction to Algorithms," 2nd ed., by Cormen, Leiserson, Rivest and Stein
# The connection matrix of the graph is represented as matrix A in the problem
# So it doesn't matter if you don't have the book. The vertexes are indexed from 0.

# The program is a free software. You may use, distribute and copy it under the 
# terms of GNU General Public License version 3.
# http://www.gnu.org/copyleft/gpl.html

# Expected output:
# $ python maxspanning.py
# vertex 0 has just been added into the growing spanning tree
# vertex 1 has just been added into the growing spanning tree
# vertex 2 has just been added into the growing spanning tree
# vertex 8 has just been added into the growing spanning tree
# vertex 5 has just been added into the growing spanning tree
# vertex 6 has just been added into the growing spanning tree
# vertex 7 has just been added into the growing spanning tree
# vertex 3 has just been added into the growing spanning tree
# vertex 4 has just been added into the growing spanning tree

from numpy import *

def spanningTree(n,A):
# neighbors[n]: // neighbors[i] is the list of neighbors of i in graph A
 neighbors=[];# neighbors[n], neighbors[i] is the list of neighbors of i in graph A
 for i in xrange(0,n):
  for j in xrange(0,n):
   neighbors.append([]);
   if( A[i,j] != 0):
    neighbors[i].append(j)

 treeNeighbor=[]; 
 # initially treeNeighbor[i] = -1;
 # during spanningTree(),
 #        treeNeighbor[i] == p  if selected[p] is true and  A[i,p]
 #            is the shortest distance from i to the current tree.
 #        treeNeighbor[i] == -1 if p is not a neighbor of any node of
 #            the current spanning tree.
 #        exception: 0 is always seleted first and treeNeighbor[0] is -1
 
 selected=[];#selected[i]=1: i is in spanning tree; 0: otherwise , initialized to be 0
 distance=[];#distance[i]: the shortest distance from i to tree, initialized to be +1
 
 for i in xrange(0,n):
  treeNeighbor.append(-1);
  selected.append(0);
  distance.append(100); # \infty for positive case, 1 for negative case

 first = 1; # a marker to mark, first =1 if the abitrary starting point is to be selected.
 for i in xrange(0,n):
  v = findMin(n,first,selected,treeNeighbor,distance);
  print "vertex", v, "has just been added into the growing spanning tree"

  first = 0;
  selected[v] = 1;    

#  add neighbors of v to treeNeighbor[] and update their distance to current spanning tree
  for i, w in enumerate(neighbors[v]):#each w \in neighbors[v] 
   #w= neighbors[v][i];
   if selected[w]==0 : #only update vertexes that are NOT on the tree now
    if A[v,w] < distance[w]:
     distance[w] = A[v,w];
     treeNeighbor[w] = v; #treeNeighbor, neighbors to current spanning tree
     #A[w,v]is the shortest path from w to current tree.
 return treeNeighbor;
 
def  findMin(n,first,selected,treeNeighbor,distance):
 if (first == 1):
  v = 0;
  treeNeighbor[0] = -1;
  return v ;
 else:
  min_distance=100; #\infty for postive case, 0 for negative case
  for v in xrange(0,n):
   if (selected[v]==0) and (treeNeighbor[v] != -1) and (distance[v] < min_distance):
   # the new node must be not on the current tree and be reachable to the tree in one hop and the shortest one to the current tree
    min_v = v;
    min_distance = distance[v]
  return min_v;

A=[\
[0, 4, 0, 0, 0, 0, 0, 8, 0],\
[4, 0, 8, 0, 0, 0, 0, 11, 0],\
[0, 8, 0, 7, 0, 4, 0, 0, 2],\
[0, 0, 7, 0, 9, 14, 0, 0, 0], \
[0, 0, 0, 9, 0, 10, 0, 0, 0],\
[0, 0, 4, 14, 10, 0, 2, 0, 0],\
[0, 0, 0, 0, 0, 2, 0, 1, 6],\
[8, 11, 0, 0, 0, 0, 1, 0, 7],\
[0, 0, 2, 0, 0, 0, 6, 7, 0]\
]

A=array(A)
R= spanningTree(9, A)

Python editing and debugging on Linux

by Forrest Sheng Bao http://fsbao.net

Python is my favorite programming language. There are many Python editors, debuggers and IDEs. Some good editors/IDE are geditor(comes with Gnome desktop by default), SPE (Stani's Python Editor) and Eric. There is a debugger with GUI called Winpdb. The default debugger is IDLE, developed by Python Software Foundation.

If I don't need to debug my program, I will only use geditor or even vim simply. If I need to debug, I will use Eric. In Eric, I have all the buttons I need on the toolbar and I can debug it without calling other programs. A trick of Eric is that I can use regular expression to hide some variables I don't want to see, as you can see in this snapshot.

The reason I mention Winpdb is because sometimes I am so lazy to hide variables I don't want to see. I just want to view certain varibales I am interested in. Winpdb provides me a way in its command console. I just need to type in "v variable_name" then I can see the value of a perticular variable.

No EOF(), the Python way

by Forrest Sheng Bao http://fsbao.net

Today, I was writing a program and gonna find a function to detect whether EOF, the end of file marker, was reached. But I failed. Not such function in file I/O of Python. Then I did some googling and realized that we don't need an EOF detection function. Take a look at this: http://groups.google.com/group/comp.lang.python/browse_thread/thread/ed25388487b3ac7b

If you try to read some stuff, but you get nothing, then you have reached the EOF. Stupid and simple, ah?

My program now looks like:

while True:
testline = f.readline()
if len(testline) ==0:
break # EOF
#do some other stuff

According to Python documents, I don't think I need to do this when using readline() function, because "Read until EOF using readline() and return a list containing the lines thus read". http://docs.python.org/lib/bltin-file-objects.html

Plotting in matplotlib/pylab

by Forrest Sheng Bao http://fsbao.net

I like Python, the simple but elegant syntax. I can't love any other programming languages more than Python. So I prefer to do my scientific computing in Python. I am very happy to find Numpy/SciPy package. matplotlib/pylab makes it possible for me to visualize my results. Neat!

Here is my first plotting script using matplotlib/pylab:

#This software is a free software. 
#Thus, it is licensed under GNU General Public License.
#Python implementation to read U. Bonn's data and plot it
#Forrest Bao, Feb. 26  

import sys,string

from numpy import *
from pylab import *

f1=open(sys.argv[1],'r')

data=[]
for line in f1.readlines():
 data.append(map(int, line.split()))

x = range(0,len(data));
lines = plot(x, data);

setp(lines, linewidth=0.5)
#shrink_to_aspect(lines, fig_aspect=0.5)
axis('tight')

grid(True);
xlabel('sampling points');
ylabel('quantilized steps');
title(sys.argv[1][:-4]);

savefig('eps/'+sys.argv[1][:-4]+'.eps')
savefig('png/'+sys.argv[1][:-4]+'.png')

#show(); 

Python implementation of Nussinov Folding Algorithm for RNA Secondary Structure Prediction

by Forrest Sheng Bao http://fsbao.net

#This software is a free software. 
#Thus, it is licensed under GNU General Public License.
#Python implementation to Nussinov Folding Algorithm 
#for Homework 7 of Bioinformatics class.
#Forrest Bao, Nov. 29 <http://fsbao.net> <forrest.bao aT gmail.com>

import sys,string
from numpy import *
from matplotlib import *

#read sequences
#sequences are stored in many lines
f=open(sys.argv[1], 'r')
seq=[]
for line in f.readlines():
 seq.append(string.strip(line));

def delta(l,m):
 delta=0;
 if l=='A' and m=='U':
  return 1;
 elif l=='U' and m=='A':
  return 1;
 elif l=='G' and m=='C':
  return 1;
 elif l=='C' and m=='G':
  return 1;
 else:
  return 0;

def buildDP(seq):
 L=len(seq);
 s=zeros((L,L));
 for n in xrange(1,L):
  for j in xrange(n,L):
   i=j-n;
   case1=s[i+1,j-1]+delta(seq[i],seq[j]);
   case2=s[i+1,j];
   case3=s[i,j-1];
   if i+3<=j:
    tmp=[];
    for k in xrange(i+1,j):
     tmp.append(s[i,k]+s[k+1,j]);
    case4=max(tmp);
    s[i,j]=max(case1,case2,case3,case4);
   else:
    s[i,j]=max(case1,case2,case3);
 return s;

def traceback(s,seq,i,j,pair):
 if i<j:
  if s[i,j]==s[i+1,j]:
   traceback(s,seq,i+1,j,pair);
  elif s[i,j]==s[i,j-1]:
   traceback(s,seq,i,j-1,pair);
  elif s[i,j]==s[i+1,j-1]+delta(seq[i],seq[j]):
   pair.append([i,j,str(seq[i]),str(seq[j])]);
   traceback(s,seq,i+1,j-1,pair);
  else:
   for k in xrange(i+1,j):
    if s[i,j]==s[i,k]+s[k+1,j]:
     traceback(s,seq,i,k,pair);
     traceback(s,seq,k+1,j,pair);
     break;
 return pair;

for q in xrange(0,len(seq)):
 pair=traceback(buildDP(seq[q]),seq[q],0,len(seq[q])-1,[])
  print "max # of folding pairs: ",len(pair);
 for x in xrange(0,len(pair)):
  print '%d %d %s==%s' % (pair[x][0],pair[x][1],pair[x][2],pair[x][3]);
 print "---";

Python script to solve N-Queen problem using minisat

by Forrest Sheng Bao http://fsbao.net

Only works for Linux!!!

#queen.py
#This software is a free software. 
#Thus, it is licensed under GNU General Public License.
#Python code to solve N-Queen problem
#generate DIMACS CNF file and call minisat to solve it.
#Usage: python queen.py 
#for SAT solver project of CS384, Logics for Computer Science
#Forrest Bao, Nov. 26  

import sys,string,os;

n=int(sys.argv[1]);

if (len(sys.argv)<2):
 print "Usage: python queen.py   ";
 print "Or, just: python queen.py ";
 exit();

if (len(sys.argv)<3):
 print "You don't specify output DIMACS CNF file name, use default I/O name- nqueen.cnf result";
 ofilename="nqueen.cnf";
 o2filename="result";

if (len(sys.argv)==3):
 print "Using your DIMACS CNF filename: ",sys.argv[2];
 ofilename=sys.argv[2];

if (len(sys.argv)==4):
 print "Using your result filename: ",sys.argv[3];
 o2filename=sys.argv[3];

ofile= open(ofilename,"w");
o2file= open(o2filename,"w");


str0="c"+str(n)+"-queen problem\nc By Forrest Sheng Bao http://fsbao.net\nc This software is a free software under GNU General Public License\nc assume number means grid\nc true means a queen in this grid\nc false means not queen in this grid\nc each line must has one queen\nc no two queens attach each other\nc so they cannot be true at the same time\n";
ofile.write(str0); 

#generator rule
for i in range(0,n):
 str0='';
 for j in range(1,n+1):
  number=str(j+n*i);
  str0=str0+number+" ";
 str0=str0+" 0\n";
 ofile.write(str0); 
#end of generator rule

#constraints on rows
for i in range(0,n):
 for j in range(1,n):
   number=j+n*i;
   for l in range(1,n-j+1):
    str0="-"+str(number)+" -"+str(number+l)+" 0\n";
    ofile.write(str0); 
#END of constraints on rows

#constraints on columns
for j in range(1,n+1):
 for i in range(0,n):
   number=j+n*i;
   for l in range(1,n-i):
    str0="-"+str(number)+" -"+str(number+n*l)+" 0\n";
    ofile.write(str0); 
#END of constrains on columns

#constraints on NW->SE diagonal
# part 1, upper bound triangle
for i in range(0,n-1):
 for j in range(i,n-1):
  number=j+1+n*i;
  for l in range(1,n-j):
    str0="-"+str(number)+" -"+str(number+l*(n+1))+" 0\n";
    ofile.write(str0); 

# part 2, lower bound triangle
for i in range(0,n-1):
 for j in range(0,i):
  number=j+1+n*i;
  for l in range(1,n-i):
   str0="-"+str(number)+" -"+str(number+l*(n+1))+" 0\n";
   ofile.write(str0); 
#END of constraints on NW->SE diagonal

#constraints on NE->SW diagonal
# part 1, upper bound triangle
for i in range(0,n):
 for j in range(0,n-i):
  number=j+1+n*i;
  for l in range(1,j+1):
   str0="-"+str(number)+" -"+str(number+l*(n-1))+" 0\n";
   ofile.write(str0); 

# part 2, lower bound triangle
for i in range(0,n):
 for j in range(n-i,n):
  number=j+1+n*i;
  if (number != n*n ):
   for l in range(1,n-i):
    str0="-"+str(number)+" -"+str(number+l*(n-1))+" 0\n";
    ofile.write(str0); 
#END of constraints on NE->SW diagonal

ofile.close();

exe="minisat "+ofilename+" "+o2filename;

os.system(exe);

exe2="cat "+o2filename;
os.system(exe2);