Thursday, November 23, 2017

MATLAB GUI Programming: Disabling Edit Uicontrols without Changing Their Appearance

Recently, I have been developing a MATLAB GUI for one of our offices that will allow them to solve complex optimization problems.  There are many inputs for the problem depending on the engineering details of that particular project, and a lot of quantities calculated from those entries  that need to be displayed near those inputs.

In previous projects, I just used the text uicontrol to handle this, but I always found the result aesthetically unappealing when a large amount of data needed to be summarized close to the user-entered data.  Figure 1 shows an example of this.  While it looks adequate for this particular case, I don't like the way it looks scaled up when there are dozens of user inputs.

Figure 1:  Using a text uicontrol to summarize information derived from user-entered data

I wanted to use the edit uicontrols to hold this information,  but by its nature the information in that box can be changed by the user.  Disabling the edit box grays out the box and text inside, which, in my opinion, doesn't look particularly good.  See figure 2.

Figure 2:  Examples of enabled and disabled edit uicontrols
Unfortunately, MATLAB doesn't have a way to disable the control without being grayed-out.  MATLAB's GUI capabilities are derived from Java Swing, and MATLAB only makes certain parts of the underlying Java object visible to the programmer.  Fortunately, one can access the handle to the Java object using the utility findjobj available in the Mathworks File Exchange.   Findjobj was written by Yair Altman who runs the excellent site Undocumented MATLAB, which largely deals, unsurprisingly, with undocumented features of MATLAB.

Once we have the object's Java handle, we can use its setEditable method to make the control uneditable.  An example is shown in the code below.


function edit_box
clear all; clc; close all;

f = figure();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Main entry uipanel  %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
entry_panel = uipanel('Parent', f, 'Title', 'Parameters', ...
 'Position', [0.05, 0.55, 0.5, 0.4]);

number_label = uicontrol('Parent', entry_panel, ...
 'Style', 'text', ...
 'String', 'Number of spheres', ...
 'HorizontalAlignment', 'left', ...
 'Position', [10 120, 100, 17]);
number_entry = uicontrol('Parent', entry_panel', ...
 'Style', 'edit', ...
 'Position', [120 120, 75, 20]);

radius_label = uicontrol('Parent', entry_panel, ...
 'Style', 'text', ...
 'String', 'Radii of spheres', ...
 'HorizontalAlignment', 'left', ...
 'Position', [10 90, 100, 17]);
radius_entry = uicontrol('Parent', entry_panel', ...
 'Style', 'edit', ...
 'Position', [120 90, 75, 20]);
radius_unit_label = uicontrol('Parent', entry_panel, ...
 'Style', 'text', ...
 'String', 'cm', ...
 'HorizontalAlignment', 'left', ...
 'Position', [200 90, 40, 17]);

density_label = uicontrol('Parent', entry_panel, ...
 'Style', 'text', ...
 'String', 'Density of spheres', ...
 'HorizontalAlignment', 'left', ...
 'Position', [10 60, 100, 17]);
density_entry = uicontrol('Parent', entry_panel', ...
 'Style', 'edit', ...
 'Position', [120 60, 75, 20]);
density_units_jlabel = javaObjectEDT('javax.swing.JLabel', ...
 '<HTML>g &#47; cm<SUP>3</SUP></HTML>');
[density_units_label, hcontainer] = javacomponent(density_units_jlabel, ...
 [200 60, 40, 22], entry_panel);

%  Create a font for the JLabel that matches the default MATLAB font for text
%  uicontrols
label_font = java.awt.Font('MS Sans Serif', java.awt.Font.PLAIN, 11);
density_units_label.setFont(label_font);

update_button = uicontrol('Parent', entry_panel, ...
 'String', 'Update', ...
 'Callback', @update_fields, ...
 'Position', [10 15 75 30]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Main entry uipanel  %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
summary_panel = uipanel('Parent', f, 'Title', 'Summary', ...
 'Position', [0.05, 0.05, 0.5, 0.4]);

volume_per_sphere_label = uicontrol('Parent', summary_panel, ...
 'Style', 'text', ...
 'String', 'Volume per sphere', ...
 'HorizontalAlignment', 'left', ...
 'Position', [10 120, 100, 17]);
volume_per_sphere_entry = uicontrol('Parent', summary_panel', ...
 'Style', 'edit', ...
 'Position', [120 120, 75, 20]);
volume_per_sphere_units_jlabel= javaObjectEDT('javax.swing.JLabel', ...
 '<HTML>cm<SUP>3</SUP></HTML>');
[volume_per_sphere_units_label, hcontainer] = javacomponent( ...
 volume_per_sphere_units_jlabel, ...
 [200 120, 40, 22], summary_panel);
volume_per_sphere_units_jlabel.setFont(label_font);

total_volume_label = uicontrol('Parent', summary_panel, ...
 'Style', 'text', ...
 'String', 'Total volume', ...
 'HorizontalAlignment', 'left', ...
 'Position', [10 90, 100, 17]);
total_volume_entry = uicontrol('Parent', summary_panel', ...
 'Style', 'edit', ...
 'Position', [120 90, 75, 20]);
total_volume_units_jlabel = javaObjectEDT('javax.swing.JLabel', ...
 '<HTML>cm<SUP>3</SUP></HTML>');
[total_volume_units_label, hcontainer] = javacomponent( ...
 total_volume_units_jlabel, ...
 [200 90, 40, 22], summary_panel);
total_volume_units_jlabel.setFont(label_font);

total_mass_label = uicontrol('Parent', summary_panel, ...
 'Style', 'text', ...
 'String', 'Total mass', ...
 'HorizontalAlignment', 'left', ...
 'Position', [10 60, 100, 17]);
total_mass_entry = uicontrol('Parent', summary_panel', ...
 'Style', 'edit', ...
 'Position', [120 60, 75, 20]);
total_mass_unit_label = uicontrol('Parent', summary_panel, ...
 'Style', 'text', ...
 'String', 'g', ...
 'HorizontalAlignment', 'left', ...
 'Position', [200 60, 40, 17]);

%  Make the edit entries in this section uneditable by accessing their Java
%  handle and setting editable to false.
jvolume_per_sphere_entry = findjobj(volume_per_sphere_entry);
jvolume_per_sphere_entry.setEditable(false);

jtotal_volume_entry = findjobj(total_volume_entry);
jtotal_volume_entry.setEditable(false);

jtotal_mass_entry = findjobj(total_mass_entry);
jtotal_mass_entry.setEditable(false);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Calback function  %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function update_fields(source, events)
  %  Read inputs from edit boxes
  N = str2num( number_entry.String );
  r = str2num( radius_entry.String );
  rho = str2num( density_entry.String );

  %  Calculate the needed values
  V_sphere = 4/3 * pi * r^3;
  V_total = N * V_sphere;
  M = rho * V_total;

  %  Update the edit boxes
  volume_per_sphere_entry.String = num2str( V_sphere );
  total_volume_entry.String = num2str( V_total );
  total_mass_entry.String = num2str( M );
end

end

This code produces the window shown if Fig. 3.  The edit uicontrols in the Summary pane are not changeable by the user, though one can highlight the values and copy them to the clipboard if need be.

Figure 3.  Example of window where several edit uicontrols are set to be non-editable while not being grayed out.
Another annoying quirk of MATLAB GUIs is the inability to customize the text in text uicontrols.  I am not sure what classes the underlying Java object derives its properties from, but whatever it is, it is very limited.  In cases where we need to format the text-- in the above example we have units raised to a power-- we can just substitute a Java JLabel in place of the text uicontrol.  MATLAB and Java default to different fonts, so I changed the font of the JLabel to match that of the MATLAB control by creating a Java font object with the appropriate settings and passing it to the setFont method of the JLabel.  The first instance of this in the above code occurs at the lines:


label_font = java.awt.Font('MS Sans Serif', java.awt.Font.PLAIN, 11);
density_units_label.setFont(label_font);

This is an example of using Java directly in MATLAB which is more a forte of the Undocumented MATLAB site so I won't dwell on it here.

Saturday, November 4, 2017

Simulating physical systems - Do not blindly trust computer output

In putting together some code to bring myself up to speed on a software library, I ran into an issue that I've seen many times over the years. You can write down the proper equations that describe the physics of your system, have an algorithm capable of solving the math involved, code it correctly, and unless you're very careful and know what's going on under the hood, still get the wrong answer.

My initial intent was to write a series of posts detailing some of the functionality of the Sundials package released by Lawrence Livermore National Laboratory (https://computation.llnl.gov/projects/sundials).  Sundials is a package for numerically solving ordinary differential equations and differential algebraic equations.  It also comes with a library for the solution of nonlinear algebraic equations.  We will go into detail on this at a later date.

I've used Sundials once or twice, but some current projects I'm working on require the speed and horsepower of C++ as opposed to MATLAB or Python.  I was making more of effort to familiarize myself with the full capabilities of the package.  To that end I was working through a couple of simple problems: a projectile motion problem with air resistance, and a damped driven harmonic oscillator.  The driven oscillator problem is nonlinear and is often used as an introduction to chaos in courses on computational physics.  Small changes in initial conditions or other parameters can lead to vastly different behaviors.

I started on the oscillator problem and grabbed an old book from my undergraduate days, Computational Physics by Nicholas J. Giordano.  I noticed something when reworking the example presented by the author.  As mentioned, the problem is a driven pendulum consisting of a mass hanging from a rigid rod of length $L$.  The angle between the rod and the vertical is denoted by $\theta$.  This is shown below in Fig. 1.

Figure 1:  Schematic of pendulum

The author considers the case where the pendulum is driven by a  periodic force with amplitude $F_d$ and frequency $\Omega$,
$$F_\mbox{drive} = F_d \sin (\Omega_d t).$$
We will also include a damping force proportional to the angular velocity,
$$F_{\mbox{damping}} = -q \omega = -q \frac{d\theta}{dt},$$
where $q$ is the constant of proportionality and $\omega$ is the pendulum's angular velocity.

The equation of motion of the pendulum is given by,
\begin{equation}
\frac{d^2 \theta}{dt^2} = -\frac{g}{L}\sin{(\theta)} - q \frac{d\theta}{dt} + F_d\sin{(\Omega_D t)}
\label{Main_de}
\end{equation}
subject to the initial conditions for our starting angle and angular velocity, $\theta(0) = \theta_0$ and $d\theta(0)/dt = \omega_0$.

We can rewrite Eq. \ref{Main_de} as two coupled first-order differential equations,
$$\frac{d \omega}{dt} = -\frac{g}{L}\sin{(\theta)} - q \frac{d\theta}{dt} + F_d\sin{(\Omega_D t)} $$ and
$$\frac{d\theta}{dt} = \omega.$$

The author solves these equations using an Euler-Cromer algorithm where we step a small time interval  $\Delta t$, and calculate approximate solutions at time  $t_{i+1} = t_i + \Delta t$ as follows,
$$\omega_{i+1} = \omega_i - \frac{g}{L}\sin(\theta_i)\Delta t - q\omega_i \Delta t + F_d\sin(\Omega_d t) \Delta t$$
$$\theta_{i+1} = \theta_i + \omega_{i+1}\Delta t.$$

Below is C++ code (a translation of the code in the textbook) implementing the above algorithm for the following model parameters.


  • g = 9.8 $\mbox{m}/s^2$
  • L = 9.8 m
  • $F_d = 1.2 \;\mbox{s}^{-2}$
  • $q_d = 0.5 \;\mbox{s}^{-1}$
  • $\Omega_d = 2/3 \;\mbox{s}^{-1}$
  • $\theta(0) = 0.2$
  • $\omega(0) = 0\;\mbox{s}^{-1}$



#include <iostream>
#include <math.h>

using namespace std;

int main()
{
        double g       = 9.8;
        double L       = 9.8;
        double F_d     = 1.2;
        double q_d     = 0.5;
        double Omega_d = 2.0/3.0;

        double period = 2 * M_PI * sqrt(g/L);
        int i = 0;
        double t = 0;
        double dt = 0.04;

        double theta = 0.2;
        double omega = 0.0;

        while (t <= 10 * period)
        {
                i = i + 1;
                t = t + dt;
                omega = omega - g/L * sin(theta) * dt - q_d * omega * dt +
                        F_d * sin( Omega_d * t) * dt;
                theta = theta + omega * dt;
                cout << t << "\t" << "\t" << theta << "\t" << omega << endl;
        }

        return 0;
}

Fig. 2 shows the results of this integration.
Figure 2.  Results of an Euler-Cromer integration of the pendulum's equations of motion with a time step of 0.04 seconds.


I mentioned this came about in an attempt to write an article detailing some of the functionality of Sundials.  For the sake of argument, let's take the output of Sundials as correct and overlay those results on the data obtained from Euler-Cromer.
Figure 3.  Euler-Cromer integration with the results from Sundials overlay ed.  The two methods agree up to around 30 seconds, then begin to drastically diverge.


OK, there is an obvious discrepancy here.  After starting out the same, the two different simulations diverge considerably as time goes on.  Is it the algorithm itself?  Well, let's consider the case with no driving force and dissipation.  Let's also take the initial displacement $\theta$ to be small such that we can make use of the approximation $\sin(\theta) \approx \theta$.  In other words, we want to solve the equation,
\begin{equation}
\frac{d^2 \theta}{dt^2} = -\frac{g}{L}\theta.
\label{simple_de}
\end{equation}
This does have an analytical solution,
$$\theta(t) = \theta_0 \sin\left(\sqrt{\frac{g}{L}} t + \frac{\pi}{2}\right).$$
Let's rerun our code with the driving force and dissipation set to zero and plot Eq. \ref{simple_de} on top of it.
Figure 4:  Comparison between analytical results fro a simple harmonic oscillator and the numerical solution using the Euler-Cromer algorithm.

Well, OK.  Those two curves are identical so the issue probably isn't with the code itself.  The issue is the time step used in the Euler-Cromer code.  As I mentioned earlier, the pendulum problem is often used as an introduction to chaotic behavior-- slight differences in parameters or initial conditions can lead to radically different behavior.  We can see that by running the sundials code for slightly differnt initial starting angles.  In the plot below, we start with our original initial position and solve for three cases where we increment the initial starting angle by 2%.  This works out to be slightly less than 0.25 degrees each time.  That isn't a big change, so we shouldn't see much of a change.  Indeed, we don't see that big of a difference between the plots until roughly 35 seconds as the four curves begin to diverge considerably.
Figure 5:  Solution to the damped driven harmonic oscillator problem for four slightly different starting angles.  The solutions are similar until around 35 seconds at which point they begin to differ significantly.

The Euler-Cromer algorithm isn't particularly accurate, so errors in the results at each time step can accumulate faster than in other, more accurate integration schemes.  This isn't really important in linear problems-- the simple harmonic oscillator plotted above is just fine-- but when we work with the nonlinear problems those small errors change the predicted results entirely.  We can confirm this is indeed what is happening by reducing the time step in the Euler-Cromer integration.  In Figure 6, we take this time step to be 0.0001 seconds.  As can be seen, this now matches well with the results out of Sundials, although some discrepancy is creeping in towards the end of the simulation.
Figure 6:  Comparison between Sundials and Euler-Cromer algorithm with the time step in Euler-Cromer set to 0.0001 S
Other issues can arise in nonlinear systems that produce erroneous results.  Changing the integration tolerance can affect the output.  Consider Figure 7 below.  These simulations were written in MATLAB using its built-in  ode45 differential equation solver routine.  Simply changing the relative tolerance the solver uses can give very different answers.  Roughly speaking, the relative tolerance sets accuracy of the solution to a given number of digits.  Figure 7 plots the solutions for relative tolerances of $10^{-3}$ and $10^{-6}$ as well as the solution out of Sundials.
Figure 7: Plot of solutions using relative tolerances of 1e-3 and 1e-6 as well as the solution as calculated by Sundials.  For clarity the density of plotted points is lower than in previous plots.

The bottom line: Even simple looking systems can produce unexpected results and be very sensitive to the algorithms (and internal settings of the algorithms) used to solve them.  Don't just blindly trust the output.  Always double check and build in various means of sanity checking your results.