Friday, February 24, 2012

Using HTML to View Large Sets of Plots - An Example in Python



This example doesn't work because of Blogger limitations. However if you run the example you will be able to select graphs in the generated HTML page.

Problem Statement

You have a program which generates lots of similar plots that end users would like to compare and explore. The end users may not be able to install any code. You can't setup a web server to nagivate the data set. You can not install any new programs on their window's desktop. How to you provide a solution?

Discussion

You can assume that any modern computer at least has a copy of Firefox, Safari, or Explorer. Since there browers support javascript (except in the worst case security settings), you can build a very lightweight data viewer using a few simple methods. The most important design decisions when generating the plots is to name the plots so they are easy to recreate from selections an end user might make.

Example

The following snippet of Python code generates 9 graphs that have random numbers plotted on two axes using different colors and markers. There are three choices of colors and three choices of markers. After generating the plots and saving them, the script creates an HTML file which simplies the navigation of the images. A user can open the HTML page and select the graph by changing the form selections at the top of the page.

There are a couple of key concepts that help make this work:
  • The plot names can be created from selections using javascript. For example, in this example there are three colors and three different markers. All of the plot file names are formed by concatenating the color and marker description to form a plot name.
  • When a user changes their choice of color or marker a javascript function rebuilds the plot file name and causes the browser to reload the image by change the img source.
  • The python script uses templates to set up the bulk of the HTML page, then substitutes in specific options for the user after the plots have been generated. 

import pylab as plt
from random import random
from string import Template

colors  = {'Red_Plot':'r',
           'Blue_Plot':'b',
           'Green_Plot':'g',
           }
markers = {'Circle_Plot':'o',
           'Square_Plot':'s',
           'Diamond_Plot':'d',
           } 

plt.figure()
for c_key in colors.keys():
    for m_key in markers.keys():
        plt.clf()
        plot_name = c_key + ',' + m_key + '.png'
        x = [random() for i in range(0,100)]
        y = [random() for i in range(0,100)]
        color = colors[c_key]
        marker = markers[m_key]
        plt.plot(x,y,color+marker,markersize=15)
        plt.savefig(plot_name)
        
HTML_template = Template('''<head>
   <script language="JavaScript"><!--
      function sel_plot() {
         // only do this if the brower supports images
         if(document.images) {
            // get plot color name
            var e=document.getElementById("color_name");
            var c_name = e.options[e.selectedIndex].text;
            // get plot marker name
            var e=document.getElementById("marker_name");
            var m_name = e.options[e.selectedIndex].text; // create the filename
            // build the filename from these selections
            var plot_filename = c_name + "," + m_name + ".png";
            // cause the correct plot to be loaded
            document["plot"].src = plot_filename;
            }
         }
      // select the plot to display initially after loading document
      window.onload = function() { sel_plot() };
      // silence errors
      window.onerr = null;
   </script>
</head>
<body>
   <center>
      <form name="Plot Select Form" id="plot_select_form">
         <select id="color_name" size="3"
            onchange="sel_plot()">
            $color_options
         </select>
         <select id="marker_name" size="3"
            onchange="sel_plot()">
            $marker_options
         </select>
      </form>
      <img name="plot" src="dummy.png" height="200"
         width="500">
   </center>
</body>
''')
        

# build the option strings
def build_options(opt_dict):
    s = ''
    for i,key in enumerate(opt_dict.keys()):
        if i==1:
            s += '<option selected>'
        else:
            s += '<option>'
        s += key
        s += '</option>\n'
    return s    

color_options = build_options(colors)
marker_options = build_options(markers)

# build the HTML page        
HTML = HTML_template.substitute(color_options=color_options,
                     marker_options=marker_options)
    

# write the HTML page
f = open('example.html','wb')
f.write(HTML)
f.close()

Test Configuration

  • PythonXY 2.7.2.1
  • IE 9


This work is licensed under a Creative Commons Attribution By license.

Saturday, December 24, 2011

How to Remove Noise from a Signal using Fourier Transforms: An Example in Python

Problem Statement:

Given a signal, which is regularly sampled over time and is “noisy”, how can the noise be reduced while minimizing the changes to the original signal.

Discussion:

A common problem in reconstructing data is elimination of noise. Noise can corrupt a signal through many means: quantization, measurement noise, errors in sampling time, sensor bias, sensor nonlinearities, signal cross coupling, etc. There are many methods which can be used to eliminate the noise on a signal. Some common approaches include use of a linear filter, Kalman filtering, Wiener filtering, construction of a custom optimization problem, and any number of ad-hoc approaches. Furthermore, the filtering of the signal can be done causally or noncausally. In a causal filter, the filtering for a particular sample if determined based only on the sample received before that sample. In acausal or noncausal filtering, all of the data before and after the sample can be used to determine the ‘best’ value for each sample.
One way to quickly filter a dataset without much effort is to use a Fourier transform. A Fourier transform is a way to decompose a signal into a sum of sine waves. The amplitude and phase associated with each sine wave is known as the spectrum of a signal. If the spectrum of the noise if away from the spectrum of the original signal, then original signal can be filtered by taking a Fourier transform, filtering the Fourier transform, then using the inverse Fourier transform to reconstruct the signal.
Lets start with a  simple example. Consider a signal consisting of a single sine wave, \(s(t)=sin(w*t)\). Let the signal be subject to white noise which is added in during measurement, \(s_{measured}(t)=s(t)+n\). Let the \(F\) be the Fourier transform of \(s\). Now by setting the value of \(F\) to zero for frequencies above and below \(w\), the noise is reduced. Let \(F_{filtered}\) be the filtered Fourier transform. Taking the inverse Fourier transform of \(F_{filtered}\) yields \(s_{filtered}(t)\). The code sample which follows illustrates these operations.

 """
Script to demonstrate the use of Fourier Transforms to acausally filter 
   a signal.

"""


__author__ = 'Ed Tate'
__email__  = 'edtategmail-dot-com'
__website__ = 'exnumerus.blogspot.com'
__license__ = 'Creative Commons Attribute By - http://creativecommons.org/licenses/by/3.0/us/''' 


from pylab import *

# setup the problem
num_samples  = 1000 # number of samples

# generate an ideal signal
f_signal  = 6   # signal frequency  in Hz
dt = 0.01 # sample timing in sec
p  = 30   # 30 degrees of phase shift
a  = 1    # signal amplitude
s = [a*sin((2*pi)*f_signal*k*dt) for k in range(0,num_samples)]
s_time = [k*dt for k in range(0,num_samples)]

# simulate measurement noise
from random import gauss
mu = 0
sigma = 2
n = [gauss(mu,sigma) for k in range(0,num_samples)]

# measured signal
s_measured = [ss+nn for ss,nn in zip(s,n)]

# take the fourier transform of the data
F = fft(s_measured)
    
# calculate the frequencies for each FFT sample
f = fftfreq(len(F),dt)  # get sample frequency in cycles/sec (Hz)

# filter the Fourier transform
def filter_rule(x,freq):
    band = 0.05
    if abs(freq)>f_signal+band or abs(freq)<f_signal-band:
        return 0
    else:
        return x
        
F_filtered = array([filter_rule(x,freq) for x,freq in zip(F,f)])

# reconstruct the filtered signal
s_filtered = ifft(F_filtered)

# get error
err = [abs(s1-s2) for s1,s2 in zip(s,s_filtered) ]

figure()
subplot(4,1,1)
plot(s_time,s)
ylabel('Original Signal')
xlabel('time [s]')

subplot(4,1,2)
plot(s_time,s_measured)
ylabel('Measured Signal')
xlabel('time [s]')

subplot(4,1,3)
semilogy(f,abs(F_filtered),'or')
semilogy(f,abs(F),'.')
legend(['Filtered Spectrum','Measured Spectrum',])
xlabel('frequency [Hz]')

subplot(4,1,4)
plot(s_time,s_filtered,'r')
plot(s_time,s,'b')
legend(['Filtered Signal','Original Signal'])
xlabel('time [s]')


show()


The output from this script shows the initial signal, followed by the measured signal which is corrupted with noise. The spectrum of this signal is filtered to recover the original signal.





References and useful links:


Test Configuration:
  • win7
  • PythonXY 2.7.2.1

Thursday, August 4, 2011

Random Notes: Validation of models vs physical systems when controls are developed independently

Problem Statement

How do you validate a model of a system against a physical system when a controller is necessary to make the system operate and the the operational policies of the controllers were developed independently.

Discussion

Consider a development process with a well defined operation cost model which drives both the model and physical system to optimization that operations cost. If the model uses full state feedback for control and the physical system uses either full state feedback or implements output feedback with state estimators and relies on certainty equivalence, there is no guarantee that under identical test conditions the trajectories of the cost and the trajectories of the states will be identical even if both system are optimal with respect to the operational costs. This is because optimal operational costs are unique, but there are no guarantees that optimal tratectories are unique. Furthermore, if the controls in both cases are not globally optimal, by only near optimal, then likelihood of non-unique trajectories is even more likely. However, because the operational costs can be unique, the validation exercise can be decomposed into two validation steps.

First, the equations which model the physics can be validated against test data on the physical system by measuring the states in the real system, then substituting the integrator in the model with the state measurements. Ideally, the physical system could execute both policies. The error in costs, and derivative calculations can be compared to quantify the error between the model of the physics and the real physics.

Second, once the errors in the model of the physics are quantified, the error in the costs under the different controllers can be quantified. Ideally, from this step, the optimality of the controls wrt to a globally optimal controller (or minimizing controller) can be established. Once interesting possibility is to use policy improvement to see if the independently developed policies can merged for better performance. Alternatively, if there are unexplained differences, then the constraints respected by the different policies need to be reconciled. Things like robustness may also contribute to differences. Robustness will in general be a driven by different view of noise and risk sensitivity. Following on that, there is the possibility that the equation structure in the different policies lead to different performance limitations. Again, policy improvement may provide a way to identify these structure imposed limitations.

This work is licensed under a Creative Commons Attribution By license.

Wednesday, August 3, 2011

Random Notes: Thoughts on metrics for engineering tools

  • If training is required for successful usage, what is the average half-life of the training. In other words, if a group of users is trained, how long until 50% of the users will forget some key aspect of tool usage which drives them to abandon the tool?
  • What is the average time for user to need to go to the help files to complete a task if they do not use the tool constantly?
  • Can a user successfully use the tool without training?
  • Are the documentation and examples sufficient for self learning?
  • How many actions are required to complete a ‘quickstart’ example?
  • How many decisions are required to complete a ‘quickstart’ example?
  • How many choices are the in each decision in a typical workflow?
  • How difficult is it to integrate the tool into automated work flows?
  • How difficult is it to customize the tool?
    • Can a power user customize the tool?
  • How long does it take to introduce a new feature in the tool?
  • How many sentences does it take to describe why a user should adopt the tool?
  • In the absence of process enforcement, would the users naturally adopt this solution?
  • What is the time saving for the individual, team, and organization from the adoption of the tool?
  • If the tool reduces error rates, is there feedback to the users to help them understand the improvement?
  • Can the input and output to the tool be reused so that the effort can be reapplied?
  • What is the ‘activation potential’ to get a new user to adopt the tool?
    • Do new users request access to the tool?
  • In a corporate setting, how difficult are the permissions to manage?
    • If a new user if not setup, will the team be able to duplicate the permissions with without calling the developers?

This work is licensed under a Creative Commons Attribution By license.

Monday, August 1, 2011

Draft: Notes on dynamic programming equations which solve cost models for dynamic systems

 

Deterministic Cost Models

Description

Cost Model

Dynamic Programming Equations

Restrictions

Finite Horizon Total Cost

\(J^{\pi}\left(x_{0}\right)=\sum_{k=0}^{K}\alpha^{k}\cdot c_{k}\left(x_{k},\pi\left(x_{k}\right)\right)\)

\( V_{k}^{\pi}\left(x\right)=c_{k}\left(x,\pi\left(x\right)\right)+\alpha\cdot V_{k+1}^{\pi}\left(f\left(x,\pi\left(x\right)\right)\right)\),\(\forall k\in\left\{ 0,\cdots,K-1\right\}  \)

\(V_{K}^{\pi}\left(x\right)=c_{K}\left(x,\pi\left(x\right)\right)\)

\(0\leq\alpha<1\)

Infinite Horizon Total Cost

\(J^{\pi}\left(x_{0}\right)=\sum_{k=0}^{\infty}\alpha^{k}\cdot c\left(x_{k},\pi\left(x_{k}\right)\right)\)

\(V^{\pi}\left(x\right)=c\left(x,\pi\left(x\right)\right)+\alpha\cdot V^{\pi}\left(f\left(x,\pi\left(x\right)\right)\right)\)

\(0\leq\alpha<1\)

Finite Horizon Shortest Path

\(J^{\pi}\left(x_{0}\right)=\sum_{k=0}^{K}\alpha^{k}\cdot c_{k}\left(x_{k},\pi\left(x_{k}\right)\right)\)

\( V_{k}^{\pi}\left(x\right)=c_{k}\left(x,\pi\left(x\right)\right)+\alpha\cdot V_{k+1}^{\pi}\left(f\left(x,\pi\left(x\right)\right)\right)\),\(\forall k\in\left\{ 0,\cdots,K-1\right\} \)

\(V_{K}^{\pi}\left(x\right)=c_{K}\left(x,\pi\left(x\right)\right)\)

\(0\leq\alpha\leq1\)

\(\left\{ x\in\chi|c\left(x,\pi\left(x\right)\right)=0\right\} \neq\left\{ \oslash\right\} \)

Infinite Horizon Shortest Path

\(J^{\pi}\left(x_{0}\right)=\sum_{k=0}^{\infty}\alpha^{k}\cdot c\left(x_{k},\pi\left(x_{k}\right)\right)\)

\(V^{\pi}\left(x\right)=c\left(x,\pi\left(x\right)\right)+\alpha\cdot V^{\pi}\left(f\left(x,\pi\left(x\right)\right)\right)\)

\(0\leq\alpha\leq1\)

\(\left\{ x\in\chi|c\left(x,\pi\left(x\right)\right)=0\right\} \neq\left\{ \oslash\right\} \)

Average Cost

\(J^{\pi}\left(x_{0}\right)=\underset{K\rightarrow\infty}{\lim}\frac{1}{K}\sum_{k=0}^{K}\alpha^{k}\cdot c\left(x_{k},\pi\left(x_{k}\right)\right)\)

\(V^{\pi}\left(x\right)+\lambda=c\left(x,\pi\left(x\right)\right)+V^{\pi}\left(f\left(x,\pi\left(x\right)\right)\right)\)

\(0\leq\alpha<1\)

\(V^{\pi}\left(x_{ref}\right)=0\) for some \(x_{ref}\in\chi\)

 

Stochastic Cost Models

Description

Cost Model

Dynamic Programming Equations

Restrictions

Finite Horizon Total Cost

\(J^{\pi}\left(x_{0}\right)=E^{W}\left[\sum_{k=0}^{K}\alpha^{k}\cdot c_{k}\left(x_{k},\pi\left(x_{k}\right),w\right)\right]\)

\(V_{k}^{\pi}\left(x\right)=E^{W}\left[c_{k}\left(x,\pi\left(x\right),w\right)+\alpha\cdot V_{k+1}^{\pi}\left(f\left(x,\pi\left(x\right),w\right)\right)\right]\)

\(V_{K}^{\pi}\left(x\right)=E^{W}\left[c_{K}\left(x,\pi\left(x\right)\right)\right]\)

\(0\leq\alpha<1\)

Infinite Horizon Total Cost

\(J^{\pi}\left(x_{0}\right)=E^{W}\left[\sum_{k=0}^{\infty}\alpha^{k}\cdot c\left(x_{k},\pi\left(x_{k}\right),w\right)\right]\)

\(V^{\pi}\left(x\right)=E^{W}\left[c\left(x,\pi\left(x\right),w\right)+\alpha\cdot V^{\pi}\left(f\left(x,\pi\left(x\right),w\right)\right)\right]\)

\(0\leq\alpha<1\)

Finite Horizon Shortest Path

\(J^{\pi}\left(x_{0}\right)=E^{W}\left[\sum_{k=0}^{K}\alpha^{k}\cdot c_{k}\left(x_{k},\pi\left(x_{k}\right),w\right)\right]\)

\(V_{k}^{\pi}\left(x\right)=E^{W}\left[c_{k}\left(x,\pi\left(x\right),w\right)+\alpha\cdot V_{k+1}^{\pi}\left(f\left(x,\pi\left(x\right),w\right)\right)\right]\)

\(V_{K}^{\pi}\left(x\right)=E^{W}\left[c_{K}\left(x,\pi\left(x\right)\right)\right]\)

\(0\leq\alpha\leq1\)

\(\left\{ x\in\chi|c\left(x,\pi\left(x\right)\right)=0\right\} \neq\left\{ \oslash\right\} \)

Infinite Horizon Shortest Path

\(J^{\pi}\left(x_{0}\right)=E^{W}\left[\sum_{k=0}^{\infty}\alpha^{k}\cdot c\left(x_{k},\pi\left(x_{k}\right),w\right)\right]\)

\(V^{\pi}\left(x\right)=E^{W}\left[c\left(x,\pi\left(x\right),w\right)+\alpha\cdot V^{\pi}\left(f\left(x,\pi\left(x\right),w\right)\right)\right]\)

\(0\leq\alpha\leq1\)

\(\left\{ x\in\chi|c\left(x,\pi\left(x\right)\right)=0\right\} \neq\left\{ \oslash\right\} \)

Average Cost

\(J^{\pi}\left(x_{0}\right)=E^{W}\left[\underset{K\rightarrow\infty}{\lim}\frac{1}{K}\sum_{k=0}^{K}\alpha^{k}\cdot c\left(x_{k},\pi\left(x_{k}\right),w\right)\right]\)

\(V^{\pi}\left(x\right)+\lambda=E\left[c\left(x,\pi\left(x\right),w\right)+V^{\pi}\left(f\left(x,\pi\left(x\right),w\right)\right)\right]\)

\(0\leq\alpha<1\)

\(V^{\pi}\left(x_{ref}\right)=0\) for some \(x_{ref}\in\chi\)

 

Risk Aware/Averse Stochastic Cost Models

Description

Cost Model

Dynamic Programming Equations

Restrictions

Certainty Equivalence with exponential utility \(J^{\pi}\left(x_{0}\right)=\underset{K\rightarrow\infty}{\limsup}\frac{1}{K}\cdot\frac{1}{\gamma}\cdot\ln\left(E^{W}\left[\exp\left(\sum_{k=0}^{K-1}c\left(x,\pi\left(x\right),w\right)\right)\right]\right)\)    
Mean-Variance      
       
       
       

 

Cost Models That don’t work or have issues

Description

Cost Model

Issues

Expected exponential disutility \(J^{\pi}\left(x_{0}\right)=\underset{K\rightarrow\infty}{\limsup}\frac{1}{K}\cdot E^{W}\left[\textrm{sgn}\left(\gamma\right)\cdot\exp\left(\gamma\cdot\sum_{k=0}^{K-1}c\left(x,\pi\left(x\right),w\right)\right)\right]\) Does not discriminate among policies
Different version of expected exponential disutility \(J^{\pi}\left(x_{0}\right)=\underset{K\rightarrow\infty}{\limsup}\frac{1}{\gamma}\cdot\log\left(E^{W}\left[\exp\left(\gamma\cdot\frac{\gamma}{K}\sum_{k=0}^{K-1}c\left(x,\pi\left(x\right),w\right)\right)\right]\right)\) Generally reduces to cost average
     
     
     

 

References

This work is licensed under a Creative Commons Attribution By license.