Tag Archives: Programming

Solar Cells and the Lambert W Function

Computations for the one-diode model for solar cells, if done using the exact formula with Lambert W function, are likely to produce arithmetic overflow or underflow. That is a constraint on the ability to implement such calculations in Fortran or C, or on microcontrollers. The solution: use a coordinate transformation of the computation problem. If the problem were being solved on graph paper, the coordinate transformation would be achieved by using log-log graph paper.

I gave a talk at a recent conference, “Celebrating 20 years of the Lambert W function”. Title: Solar Cells and the Lambert W function. Joint work with my colleague S. R. Valluri. The slides are available at Researchgate, at this URL:

https://www.researchgate.net/publication/305991463

Best wishes,
Ken R.
11-Aug-2016

Finding Good Recommendations

Leave a reply

There are many online services which attempt to recommend movies, books, webpages etc that someone will like. In some implementations, recommendations are based upon finding new items that are liked by other users: If you like A, and someone else likes A and B, then B is perhaps a good recommendation for you. The difficulty is that establishing your profile can be a tedious task, as you have to initially indicate several items that you like. On the order of twenty items, perhaps.

A new algorithm, developed by Evgeny Frolov and Ivan Oseledets of the Skolkovo Institute, provides a much less time consuming, and likely more accurate, way of establishing your preferences. It uses information about items that you do not like, as well as about items that you like. Roughly stated, if you do not like item B, and those who like B also like C, then item C is perhaps not a good suggestion for you.

The details of their algorithm are subtle, and designed for efficient operation. It is not just graph searching. See Arxiv 1607.04228 for their paper — linked below — and a press release also linked below.

Best wishes,
Ken Roberts
01-Aug-2016

http://arxiv.org/abs/1607.04228
Evgeny Frolov and Ivan Oseledets — Fifty Shades of Ratings: How to Benefit from a Negative Feedback in Top-N Recommendations Tasks

Press Release — Skoltech scientists have created an algorithm that improves the quality of recommender systems

Raspberry Pi 3 Multi-Computer

Leave a reply

It now makes sense to use Raspberry Pi 3 boards to build a multi-computer for serious scientific work. In a previous post on 25-Jan-2016 I described my tests using a Raspberry Pi 2 for density functional theory (DFT) calculations using the ELK software. The RP 2 operates at about 1-10th the speed of a Laptop purchased for scientific work, the laptop having an Intel i5-core processor and 4 GB of memory. Details in the previous post. Bottom line: RP version 2 is not economic for building a multi-computer for the ELK DFT tasks.

However, I have now had an opportunity to test the Raspberry Pi version 3, released earlier this year. It is faster, and completed the ELK calculations in 43 minutes elapsed time, vs 79 minutes on the RP version 2, and vs 8 minutes on the i5-core laptop. That means it takes 5.5 RP3 boards to equal the capability of the laptop. Since each RP3 (board only) costs $35 Cdn from a typical supplier (or $45 if one includes a power adaptor, or $75 if one wants a case, cables, etc as well), one can set up a config of about 6 RP3 boards for something like $250. The laptop costs $300 as a refurb system.

There are of course complexities. One has to mount the RP3 boards somewhere — a bit of lumber should suffice for home brew. And six power supplies can probably be replaced by a single supply of larger amperage capacity. Ethernet cables are needed, and a switch, but most tinkerers have that sort of stuff around. Still, all considered, one can probably build a very nice multi-computer out of RP3 boards at about three-quarters the cost of an equivalent multi-computer based upon i5-core refurb laptops.

I’m not planning to actually build this multi-computer. I already have what I need for my calculational tasks. But it’s nice to see the Raspberry Pi become suitable for this serious calculation.

Best wishes,
Ken Roberts
05-May-2016

ELK Software and Raspberry Pi

Leave a reply

ELK is a software package for density functional theory calculations. Raspberry Pi is an inexpensive computer. I was interested in whether the RasPi might be useful for ELK calculations. This post reports on some timing tests.

ELK is open source software, written in Fortran, available via the SourceForge link given below. ELK runs its calcs on multiple parallel processors. My typical config for using ELK is to run on a two-core or four-core processor such as an Intel i3, i5 or i7 based Laptop, with Linux (Slackware) as the operating system, and OpenMP as the task coordination mechanism. I usually run several jobs, with various parameter choices, exploring some material model. Each job might take from a few hours to a couple days, depending upon the fineness of the grid used for modelling. On a supercomputer cluster, I have in the past run up to 100 jobs concurrently, with a selection of parameter choices, in order to develop an understanding of how the parameters affect the calculation results for the material.

More recently I have been running jobs sequentially, using a more directed exploration of the parameter space vs model results. That latter approach allows for more interaction with the investigation as it progresses, and leads to better intuition. One of my beliefs is that the dynamic tension between improving calculation speed and improving model and math, leads to better understanding. That belief goes back decades, to some success I had in number theory problems by using the tension between calculation and description. To do the work faster, instead of using a faster computer, one can replace parts of the problem description with better math. Eventually, sometimes, the problem collapses into good math and a fairly rapid calculational model. I have a hope in the back of my wish list that certain DFT, QCD, and Molecular Dynamics problems will someday be discovered to have such simplifications.

Anyway, given that approach, I was interested in perhaps using the RasPi as an ELK machine, thereby leaving my main Linux laptop free for other work. The RasPi (couple of links below) has impressive specs, with the model 2B, the current product, having a quad-core processor, 1 GB memory, 8 GB to 32 GB of local disk (depending on the sD photo card installed) with net 4 GB to 28 GB of free storage for user files, an ethernet port, four usb ports for keyboard, mouse, etc, and … best of all … very modest cost. RasPi comes with Linux (Debian based) and a copy of Mathematica. Although I don’t use Mathematica at present, having “paid my dues” already by learning Maple, and not wishing to have to re-learn the quirks of a new tool unless necessary, the low cost access to Mathematica is a nice plus, of possible future benefit.

So, I ran some timing tests. The results in summary: ELK installs on RasPi with no difficulty (details below) and runs its post-install tests successfully. However, RasPi is slow in running ELK, taking 79 elapsed minutes to complete the collection of 18 test calculations. The “top” and “1” command sequence shows that all four of the processor threads were active, so the problem may lie with some other aspect — for example, speed of the sD card or the low (1 GB) memory resulting in less in-memory-buffer file access. Whatever. In contrast, my Intel i5-core based laptop (Lenovo ThinkPad T420 with disk replaced by a 1 TB drive, and Slackware installed instead of Windows) completes the test calculations in 8 minutes. With a 10x speed difference, the ELK approach is not cost effective. It might cost about $600 to get ten stripped down RasPi boards (with board itself, with only power and ethernet connected, and controlled via ssh login sessions), whereas the T420 laptop cost only $300 refurbished. So, I will stick with laptops. However, I look for the next generation of RasPi boards. I think it is an excellent product. There are other applications besides ELK. GROMACS, for molecular dynamics calculations, for instance, might be suitable for the RasPi. For the future.

Some details about the RasPi setup may be helpful to record, for others. The RasPi in Canada is available from canakit.com and one can find other suppliers elsewhere. I was pleased with the speed of delivery from Canakit. I got the full kit, eg with Wifi dongle, but one just needs a stripped down kit for second and subsequent RasPi boards. I replaced the 8 GB sD storage card with a 32 GB sD card, with the latest operating system version (Nov-2015) downloaded from the RasPi support website. The OS boots to an X-Windows screen, and one logs in as user “pi”. To get root access and a command line, I used “sudo passwd root” and set my own password on the root user. Then one can work a bit easier. The aptitude utility was used to install the software needed for ELK, including “aptitude xxx” commands where xxx was these in sequence: “update”, “search fortran”, “install gfortran”, “search openmp”, “install libblacs-openmpi1”, “search lapack”, “install liblapack3”, “search fft”, “install fftw3”. ELK was set up to use the GNU fortran compiler (gfortran), etc.

It is interesting that the ELK test runs produce some warning messages regarding IEEE floating point arithmetic errors: IEEE_INVALID_FLAG, IEEE_OVERFLOW_FLAG, IEEE_UNDERFLOW_FLAG, IEEE_DIVIDE_BY_ZERO. These exceptions are discussed in a post at SourceForge, URL given below. Although I found the specific recommendations of that post to be ineffective, the insight into cause, and why it “does not matter”, is useful — see the third comment in that post. It is something that can guide investigation. Running ELK on Slackware / Intel i5-core with gfortran does not produce those errors, though it produces identical test results. I think that is simply because floating point exception warnings have been turned off in the latter configuration.

That’s it about RasPi and ELK. Each device/program is worth one’s attention if it matches one’s objectives, as these do with my objectives. Neither is perfect, but they are understandable, affordable, open for improvement, and a net benefit for serious work.

URLs below.

Best wishes,
Ken Roberts
25-Jan-2016

For ELK density functional theory calculational software:
http://sourceforge.net/projects/elk/

For Raspbery Pi supply, and for community projects website:
https://www.element14.com/community/welcome
https://www.raspberrypi.org/

Discussion of ELK tests producing IEEE floating point exceptions:
http://sourceforge.net/p/elk/discussion/897820/thread/e87237ad/

[end]

Solar Cells 2

Leave a reply

A solar cell produces an electric current when its material surface is illuminated. In order to standardize the methods of measuring solar cell performance, and comparing different cells, certain conventions are adopted. For instance, standard illumination, which is the amount and wavelengths of the light falling on the solar cell’s surface.

A typical solar cell might have an area of 2.2 cm^2. Sometimes the current produced by a solar cell is quoted as so-many amps (or milliamps or microamps) under standard illumination. That statistic is appropriate when considering a particular solar cell or comparing two models of solar cells. Other times the current produced by a solar cell is quoted as a “current areal density”, that is, so many amps per square centimeter of solar cell material. That statistic is appropriate when considering a proposed solar cell material, or comparing two types of material or two choices of processing options for preparing solar cells using a particular material.

Solar cells are assembled into modules, which are groups of individual solar cells connected together (eg, in series) pre-packaged for easy handling. Modules are built into arrays, which make up solar panels; the panel includes the structural framework. Solar panels are the structures one notices on roofs and in fields.

The context of our recent work has been a mathematical model which is used for describing individual solar cells. As you can see there are many other aspects of solar cells and solar panels which are worth investigation (and have been investigated). Our work is just one part of the efforts of thousands.

I’m no expert on solar power. But I want to share the bits of information (or misunderstanding!) which I’ve gathered during my work. In the rest of this post I’ll describe some general aspects of solar cell performance. The details of the math model I’ll discuss in other posts — or you can go directly to the working paper at Researchgate if you want.

There are two measurements on the performance of a solar cell which are relatively easy to make: Isc = the short circuit current, and Voc = the open circuit voltage. Isc is determined by the amount of current the solar cell will push when it is under standard illumination and its contacts are shorted (no voltage difference); that current is largely determined by the series resistance within the solar cell. Voc is determined as the voltage difference between the contacts of the solar cell, again under standard illumination, when the load is missing — ie, an open circuit (no current flowing), or the voltage into an infinite load. Under other loads, there is a current less than Isc, and a voltage difference less than Voc.

The “I-V characteristic” for a solar cell is the curve which relates the current and the voltage. The graph shows a typical example, for a silicon solar cell manufactured about 3 decades ago. The current is shown as a negative number because of the convention used in the test circuit. What you can see in this curve is that this solar cell has Isc about 102 milliamps and Voc about 520 millivolts. The solar cell puts out a fairly constant current, between 102 milliamps and 95 milliamps, regardless of the load resistance, as long as the voltage drop does not exceed say 420 millivolts. That is, the load resistance can be up to say 420/95 = 4.3 ohms. As the load resistance goes higher, the solar cell’s current output drops rapidly towards zero. A load resistance of about 4.3 ohms, for this solar cell, is the “sweet spot” when the I*V product is a maximum. Power equals I*V (for direct current circuits like a solar cell), and hence that sweet spot is known as the “maximum power point” or MPP of the solar cell. The “fill factor” or FF of the solar cell is the ratio between the I*V product at the MPP, and the Isc*Voc product.

The maximum power point can shift if the illumination on the solar cell changes. For example, if the solar cell falls into shadow. The calculations to determine the MPP for an array of solar cells, for load balancing, can represent a lot of work. The MPP can be estimated by searching, varying the load slightly, but sometimes the power curve will have multiple peaks, and automated searching may not find the optimal operating conditions. It is important to have a good mathematical model for the I-V characteristic, in order to rapidly find the MPP under changed illumination. Further, that model should be suitable for implementation on very inexpensive microprocessors, such as may be used in a field installation. That’s where the method which S. R. Valluri and I recently described in our working paper is likely to be most useful.

That’s it for this post. More later!

Best wishes,
Ken R.
20-Dec-2015

BB Playbook as Gutenberg Book Reader

Leave a reply

I wanted to make a lightweight e-book reader with minimal expense. An old Blackberry Playbook sitting around, plus free Project Gutenberg e-books, turned out to be a zero cost method. I’m going to describe the details here, because they may be useful to someone else also.

The Playbook has a nice app “100 Free Ebooks” which has the top 100-some Gutenberg ebooks downloaded onto the Playbook’s ramdisk, ie there is no need for net access when reading. That is ideal for travel, because there are no expiry dates on the books, and the internet is not used while on the road. However, how to add more of the Gutenberg books? There are almost 50,000 of them, so the top-100 most frequently downloaded is only 0.2 percent of the whole library. I wanted some of the books which are personally interesting, though not in the top 100.

The solution is this:

1) Install Blackberry Desktop (BBD) software on a Windows machine (Win-7 in my situation), and connect the Playbook (PB) via USB cable to the Win7 system. Use BBD or direct disk access via Win7 to create /books/gut directory to receive the additional Gutenburg book files.

2) Download html (or sometimes, pdf) files for the Gutenberg books from their website http://www.gutenberg.org. The files can be saved directly into the Device/books/gut directory on the PB.

3) On the PB, use browser to go to file://accounts/1000/shared/books/gut/ which is the pathname where the books were stored, in terms of the browser’s view of the file system. Put a bookmark on the home screen, and call it something like “Gut Books”.

Here’s a picture of the resulting screen. One just picks the desired book from the list, and reads it as an html or pdf file, using either the browser or Adobe reader respectively.

It would also be possible, presumably, to prepare an index.html file which would give a nice-looking front end to one’s personal ebook collection. That’s basically what the top-100 application does.

Best wishes,

Ken Roberts

11-Aug-2015

Disable Apps to Improve Performance

Leave a reply

For several months now, I’ve been doing some fairly heavy duty calculations, finding roots of polynomials. I got a new laptop, just for that purpose — it sits at home and crunches numbers.

One of the amusing things about it, is that it has a lazy operating system (Windows 8.1). As the laptop crunches along, running three copies of the calculation software in parallel, from time to time the operating system offers helpfully that it can disable some apps (= applications) to improve performance.

It’s a classic instance of narrow-mindedness on the part of the operating system development team. For them, performance is measured by the amount of unused capacity in the system. it apparently does not occur to them, at a top-of-mind operational level anyway, that performance for the user (me) might consist in getting useful work completed.

Another recent activity has been to watch some episodes of Mr. Selfridge, about the founder of the London department store, Harry Gordon Selfridge, and to read the predecessor book upon which the TV soap is based. What made HGS so successful was an unwavering commitment to service to the customer. He would figure out what the customer wanted (even if the customer did not herself or himself know, ie articulate that need), and provide it — for instance, a place to have tea in the store, or to relax, etc. HGS was the originator of the slogan “The customer is always right.”

I doubt that HGS would consider that his customers should disable their activities and commitments in order to “improve performance” in their shopping!

Best wishes,

Ken R.

24-Oct-2014

Conformal Mapping Website

Leave a reply

There’s a very nice conformal mapping book and collection of mappings online at
http://harald-dalichau.de — provided by Prof. Harald Dalichau. I’m finding it a useful resource for some conformal mapping explorations I’m presently doing. He has written a book, the complete text of which (in English or in German) is online. He also provides a comprehensive dictionary of mappings, with illustrations, and also Basic programs which do the relevant calculations. It’s quite helpful.

To get the book, etc, go to the home page above, click on Bucher, then click on one of the “Conf Mapping – Chapters” links on the left of the page, or the “Conf Mapping – Mappings” link for the mappings themselves.

Enjoy!

Best wishes,
Ken Roberts
06-Jun-2014

Transcendental Equation Solution

Leave a reply

The equation $a\,e^{bx} + bx = c$ looks hard to solve. It’s said to be transcendental, ie not algebraic, like a polynomial equation. But the Lambert W function is the right tool for solving that equation exactly. Many transcendental equations are not hard to solve. This one is easy, if you have the right tools for the purpose.

The equation a*exp(bx) + bx = c above shows up in models of diodes or photovoltaic cells, in the relationship between current I and voltage V. If you solve for either I in terms of V, or for V in terms of I, you will end up with an equation of the above form. See the website http://pvpmc.org for some background about photocells and diodes, a circuit diagram for a diode equivalent of a photocell, and the Shockley model of the diode. Or there are many other resources.

You might think that equation above is a bit special, because the constant b appears in two places. But if one has some other coefficient in front of x, simply multiply the equation by a constant factor so you get bx both in the exponential and as the linear x term. And put the constant c on the right hand side. That is a very tidy way to write the equation relating a linear variable and its exponential, and that is the form which I will solve here.

Write y = exp(bx). Thus bx is log(y). Rewrite the equation as ay + log(y) = c.

We want log(ay) on the left, to match the ay term. So add log(a) to each side:

ay + log(y) + log(a) = c + log(a)

Since the sum of two logs is the log of the product, that becomes

ay + log(ay) = c + log(a)

Now take exponentials. We get

exp(ay)*ay = exp(c + log(a)) = a*exp(c)

The left hand side is in the form that defines the Lambert W function, namely w*exp(w) = z, which is the definition of w = LamW(z). So we have

ay = LamW(a*exp(c))

Divide by a, and recall that y is exp(bx), to get

exp(bx) = y = (1/a)*LamW(a*exp(c))

Now take logs, and divide by b, to get the solution for x:

x = (1/b)*log((1/a)*LamW(a*exp(c))

This solution, in its diode application, is used in modeling of diode and photocell arrays. It is very fast compared to having to find solutions by manual searching, because computer software is widely available to calculate the LamW function. (Just as there are sine, tangent, exponential, log, square root, and many other functions available in the software that we can use without having to write our own subroutines.)

Now, let’s think a bit more widely. If two diodes are in parallel, their currents add. That is like adding two Lambert W functions. If two diodes are in series, their currents satisfy a harmonic mean relationship, of the form 1/x1 + 1/x2 = 1/x. Will it surprise you that there is a LamW interpretation of that? Basically, one can make an analogue computer out of diodes, that will do arithmetic using LamW functions as a basis.

Recall the delay differential equation solution mentioned in a prior post? Also LamW. Actually, a sum of LamW evaluations, one per branch, if one is working in the complex plane. Like a Fourier series, but better because it incorporates delay. Delay of course happens in electronic circuits — more LamW interaction with analog computer ideas.

A tremendous amount of money is spent running particular programs on digital computers. Runs of models in density functional theory (DFT) or lattice quantum chromodynamics (QCD) can take days or weeks of supercomputer time. And the job is never done. Once a solution is obtained, one decides to use more points in the grid or lattice, and to use more digits of precision to deal with noise in the result arising from discrete approximations to nature.

Meanwhile, nature manages to compute quite well physically. Want to know the lattice constant for the Gold FCC crystal? Measure a chunk of the metal, determine its specific mass, and do a little arithmetic. Or use some sort of diffraction technique. Lots of ways to get the job done. Want the most cumbersome, slowest method? Use DFT. You will be lucky to choose an algorithm and discretization that gives you a value within 3 percent of reality. It is acceptable in some circles to get a rough idea, and even to use that rough idea to confirm a hypothesis about crystal structure. Maybe not for Gold — that is pretty well known — but for more complicated materials.

Nature, however, is not so tolerant. If you take a chunk of gold metal, and increase its lattice constant by 2 percent, it will melt! To get the lattice constant up by a couple of percent, you will have to add energy, and bring the metal up to about 1350-1400 Kelvin. Whereupon it melts. It’s not melting because of the temperature, not directly. Rather, at a lattice constant about 1.02 times room temperature lattice constant, gold can no longer hold its FCC crystal structure, and the atoms start to become mobile — what we call melting or liquid.

If we cannot answer all questions about nature using digital computers, maybe we should be giving nature a chance to answer some of those questions for us. Thus my fondness for analog computation, geometric solutions and such.

Best wishes,
Ken Roberts
02-Jun-2014

Vepstas Regions 4

Leave a reply

Further to the topic of Vepstas regions, defined for positive R as V(R) = the set of complex z for which f(z)=z^2/(z-1) has magnitude less than R. Linas Vepstas, in his software package ANANT for arbitrary precision number theoretical calculations, implements a method of calculating the polylogarithm of complex order s of complex argument z. He defines a complex data type cpx which is simply a pair of pointers to two arbitrary precision real values, ie the real and imaginary parts of a complex number. The routine to calculate polylog(s,z) is cpx_polylog in the mp_polylog.c module of ANANT.

Looking at the code for that module, one sees that, after some preliminary qualification of the argument, the evaluation for z in the unit circle (or on its boundary) comes down to a very neat algorithm. Here, to start things off, is a picture of the boundary of the V(1.225) region (1.225 is approximately square root of 1.5).

If z is in V(1.225), then polylog(s,z) is evaluated directly using an algorithm due to Peter Borwein, which I’ll call the PB algorithm. Otherwise, if z is in the unit circle U but not in V(1.225), a “duplication formula” for polylog is used to transform the calculation of polylog(s,z) to that of a linear combination of polylog(s,-z) and polylog(s,z^2). Suppose that z is complex, in U-V(1.225). Then z has positive real part at least 0.6, and so -z has negative real part, so lies in V(1.225) and polylog(s,-z) is readily calculated with the PB algorithm. Moreover, z^2 is at twice the phase angle of z, so swings up (or down) from the real axis. If z^2 now lies in V(1.225), then the PB algorithm can be used to calculate polylog(s,z^2), and the duplication formula can be used to obtain polylog(s,z) as a linear combination.

Otherwise, suppose that z^2 lies in U but not in V(1.225). Then use the duplication formula again, to express polylog(s,z^2) as a linear combination of polylogs of -z^2 and z^4. The phase angle of z^4 will be twice that of z^2. And so on. Eventually, one gets z^(2^n) which lies in V(1.225), and can complete the calculation. The whole process is a very nice illustration of how recursive coding can simplify a programming task conceptually.

Perhaps in another post I’ll discuss how polylog(s,z) is calculated for z outside the unit circle.

I guess I’d better say something here about the situation of real z between about 0.6 and 1.0. For z less than 1.0, successive applications of the duplication formula will produce z^2, z^4, etc tending towards the origin, so eventually get some z^(2^n) within V(1.225). For z equal to 1, the function polylog(s,1) is identical to the Riemann zeta function, so one can use some other code to calculate Riemann zeta.

Some references so you can follow up on this yourself:

Vepstas paper: Available at the URl http://arxiv.org/abs/math/0702243

ANANT code: Can be downloaded from https://launchpad.net/anant

Best wishes,
Ken Roberts
21-Apr-2014

Random vs Pseudo-Random

Leave a reply

When doing computational experiments, such as Monte Carlo simulations, one often uses a random number generator provided by some software package. That is ok most of the time, but sometimes it leads to strange correlations. I would like to tell you about one such situation. An investigation related to that situation was my first project after graduation from university, in 1968, and the origins of the problem are still fresh in memory as a cautionary tale.

At that time, IBM were the dominant supplier of mainframe computers for scientific work, and they provided a “Scientific Subroutine Package” (SSP) as a standard library of Fortran-callable procedures for matrix arithmetic and other calculations.

The SSP library included a random number generator, which functioned in a manner which was standard at that time: Given a starting integer number (seed) S, and an integer multiplier M provided by the subroutine code, each call to the random routine calculated M*S, discarding some low order bits, and some high order bits, and storing the result back into S. That result was normalized to be a floating point number in between 0 and 1, and was returned to the calling program as a “random” value. Such algorithms are called “pseudo-random” because they are deterministic, and their output only appears random if one does not know the algorithm.

An article appeared in the scientific literature describing a problem with the IBM SSP random routine. The authors had been doing population simulations, with two independent traits each of which appeared in 2 percent of the population. They would select for individuals for which two successive random variables X1 and X2 were each less than 0.02, and then calculate two more random variables X3 and X4 to simulate the presence of two other traits in those individuals. The problem was that, using IBM’s SSP routine, if X1 and X2 were each less than 0.02, that meant there were a block of zeros in the middle of two successive S values, and that in turn resulted in a very high correlation between the values of X3 and X4. A scatterplot of X3 vs X4, selecting only for cases when X1 and X2 had been less than 0.02, showed the X3 and X4 values occupying a more or less diagonal strip across the unit X3 by X4 square, an area which was only perhaps 1/10-th of the whole square. The resulting population simulation produced some decidedly unusual results!

That SSP critique article was published, but I’ve long lost the reference – sometime in 1966 or 1967, I suppose. My first job, in 1968, was to determine whether our university’s random number routine (not based upon IBM’s SSP, but instead on an algorithm from another university) had, or did not have, a similar bias when computing four random numbers successively. And, more generally, to verify that our pseudo-random generator did not have a variety of other properties that might have been a concern if used for other work. Our routine turned out to be fine — that is, it did not exhibit the patterns and problems seen within IBM’s SSP random routine. Of course any routine based upon the “multiply integers keep only middle bits” technique has some patterns, and there are some interesting aspects, both number theoretic and practical, in developing and testing a random routine. That study was called “Investigation of UC Random”, published about 1968 as a report of the University of Chicago Computation Center. It is outdated; there are nowdays much more sophisticated ways of validating pseudo-random number generators.

Returning to the question of random numbers provided by software packages. How “good” are they? Do you know? In quick exploratory work, I generally just assume they are “good enough” and carry on, focused on the outputs. Sometimes there are surprises. In a well-known, expensive and widely used package, I recently noticed that the package’s call, which was supposed to cause it to start its random number sequence at a different seed, was not effective. That is, one would get the same “random” sequence on every run, despite having written the recommended code to randomize the starting point. Costly software is not an assurance of quality software! One must be careful. Experimentation and theoretical work should be done with humility, and humility is not the predominant characteristic of successful marketing of software. The original quality work done in early versions may have been coded around, or rendered ineffective, by generations of maintenance.

Bottom line: Try your experiments several times, with different equipment, or different random number generators, or at least different starting seeds. Ideally, get someone else to write an independent program, or whatever is needed to be reasonably confident that what you are observing is derived from the phenomenon you believe you are modelling, and is not just an artifact of your model.

Randomly yours,

Ken Roberts
05-Feb-2014

Mistakes — How to Learn from Them

Leave a reply

Making mistakes is one of the best ways to learn. Yet, we are taught to believe that making mistakes is wrong, and students are penalized for mistakes. That is unfortunate, and in this post I wish to argue for the merits of making mistakes.

Watch a baby learn — moving about, grasping, throwing, speaking, and so forth. A parent would be justified in serious concern if their child did not make constant, though unsuccessful, attempts to move about, grasp, throw, or speak. The child learns to succeed through repeated attempts, most of which are unsuccessful. Thousands of attempts to move about, to grasp, to throw, or to speak strengthen the neural and motor structures which allow for later success.

Youth is the period of our lives when we learn the fastest. I think that is partly related to it being the time when we attempt the most, that is, make lots of mistakes. Perfection and predictability are in opposition to the ability to learn.

I heard a story of an engineer in a research institute, who, when he had an idea, would try it out on his colleagues, one by one. He would walk into someone’s office, set forth the idea, and it would be shot down, full of flaws. He would revise his idea, and go to the next colleague — same result. After he had visited many colleagues, his idea would have become refined by the various circumstances and considerations set forth by his colleagues. If his idea survived, it would be robust. Most ideas would not survive unscathed; in fact, most might not survive the chain of validations at all. But the good ideas — the few percent that were really worth pursuing — would become much better. That engineer must have been known to his colleagues as a person who made lots of mistakes. It takes guts, or humility, to appear “silly” in order to keep learning new material, since that invariably involves mistakes.

I was trained as a mathematician, a field which, more than most others, puts a high premium on error-free problem solving, the avoidance of mistakes. When I became a computer programmer, it took a while for me to discover that the most effective way to program is to develop partial solutions, and then refine them (debug them) against realistic requirements. Effective programming involves the art of making mistakes, repeatedly, and learning from and extending the software based upon the experience of those mistakes.

Having returned to the academic world, in retirement, as a student, I notice that academic life puts a high value on being relatively mistake-free. That is a trap, a social pressure that can lead to reduced productivity. Mistakes, lots and lots of mistakes, whether in academic or industrial settings, are important in order to make progress.

Here’s hoping that you are able to enjoy making mistakes, as much as I do!

Ken Roberts
28-Jan-2014

Exoplanet Detection and Gaia

Leave a reply

Gaia will help to detect many more exoplanets. There are many challenges — the more detailed the data, the more likely that it will be necessary to analyze planetary systems with multiple planets. Those are of course the sorts we are most curious about.

Two papers which discuss the probable effectiveness of exoplanet detection using Gaia’s data are these:

A. Sozzetti, S. Casertano, et al, 2007, “Testing Planet Formation Models with Gaia MicroArcsecond Astronomy”, http://arxiv.org/abs/0711.4903
This is a fairly early paper, good for illustrating the state of the art during Gaia planning. Awareness of that paper came via the Gaia DPAC newsletter number 04 (April 2009) — I previously gave a link for those newsletters.

A. Sozzetti, 2012, “Astrometry and Exoplanet Characterization: Gaia and Its Pandora’s Box”, http://arxiv.org/abs/1201.3484
This paper gives a more recent view of the situation including a discussion of the various software algorithms which are to be used for exoplanet orbit determination.

There are many other papers, prior and since, and we can expect an abundance of papers once the Gaia observations start to flow forth in volume. The detection algorithms will improve as they are tested against reality. The old observational data will be reprocessed, again and again. One really can do experimental science by just sitting and thinking about what one has already seen!

One interesting aspect of the Gaia project will be to trace the improvement in techniques over the life of the project, which has been in the works for over a decade. There is guidance for other long term development efforts of all sorts.

In my younger days, we had a saying, “Keep on Trucking”! I would like to suggest this for Gaia reprocessing:

Keep on Crunching !

Ken Roberts
26-Jan-2014

Gaia the Newborn Astronomer

Leave a reply

The original proposal for Gaia’s data processing and analysis consortium, from 2008, has a good explanation of the scientific aspects. It is a 700+ page document, of which the first 240-ish pages are the most interesting to read. The document can be accessed via the web page at URL http://www.cosmos.esa.int/web/gaia/dpac/consortium

Gaia’s logical structures for data processing and analysis remind me of an infant’s learning process. Gaia must be self-calibrating, as no prior star catalogue has as much precision in measurements as Gaia will produce. Later observations by Gaia will refine earlier observations and object recognitions. Gaia’s detectors have an object detection capability which ignores dark patches — comparable to edge detection in the human eye. I expect that one of the side benefits of Gaia will be significant advances in artificial intelligence.

Welcome to the newborn Gaia! She is destined to become one of earth’s foremost astronomers, a pioneer for her species.

Ken Roberts
20-Jan-2014

Gaia DPAC Newsletters

Leave a reply

DPAC is the Data Processing and Analysis Consortium for the Gaia spacecraft mission which will conduct a new star survey. There is a wealth of detail in the DPAC newsletters. The newsletters are intended both for the many people involved in the Gaia project, and for other interested people who like a bit of substance beyond the press releases.

The newsletters can be obtained at the URL http://www.cosmos.esa.int/web/gaia/dpac/newsletter and there have been 22 issues so far, over the past 6 years. I’m presently enjoying reading back issues. Start with number 1 to get introduced to the technology and terminology used in the Gaia project.

The vision within the Gaia project, and the detailed project planning including techniques for future data access, are really impressive.

Ken Roberts
20-Jan-2014

MyCloud Experiences

Leave a reply

Recently I’ve been experimenting with a WD MyCloud network attached storage (NAS) device. It’s involved a bit of exploring around and here are some of the things I’ve found out that may be of use to some other techie. I must emphasize that this is warranty-voiding modification of the NAS device.

Observations:

— MyCloud now does an analysis of files on the /shares area, to categorize them as media, and maybe a content check too judging by the time taken. I was, on a previous MyBookLive NAS, using the shares area to backup up large groups of data files. No more, as the analysis overhead trashes the NAS performance. Instead, I set up a distinct top-level directory on the DataVolume partition, and back up to there. Also it was desirable to modify varous config files to exclude utilities like wdnotify and wdmcserver from processing the personal content area.

— Suppose one turns on, then turns off, then turns on a particular category of file, let’s say Documents, in order to focus the backup on Music, for instance. The result is that previously backed up documents are flagged as “deleted” — you can see that in the recover display for those files, as they will have a horizontal line drawn thru their names. The files are still present on the backup device. But, turning Documents backup back on, and doing a backup, does not clear the “deleted” flag.

— It’s still not clear to me that the WD software (SmartWare) is properly backing up modified files. In fact, I would say it is clear to me that there is a problem with triggering backups if a file is modified. My reaction to that, aside from curiosity, is to ensure that there is another mechanism (rsync to a distinct backup area, done outside of the SmartWare utility) to ensure that there is a backup of a modified file.

— The SmartWare backup performance is slow, far slower than rsync. That is ok if it is reliable. But not at present in my opinion.

— However, on the plus side: the box itself is very nice. I think it can be used as a storage device. Cost is moderate, about $60/terabyte for the 3-TB version at list price, and it has decent performance, in fact very good performance, as an rsync puller. Warranty-voiding as mentioned, but it is quite economical compared to other makeshifts for producing a portable backup unit (such as putting a large disk drive into a discarded laptop).

— Maybe I should mention security considerations briefly. I will not say much explicitly. Basically, not good. It is advisable to devise some techniques to protect your data, not rely upon the config of the MyCloud.

Hope that is helpful, and good luck!

Ken Roberts

10-Dec-2013

Blackberry Playbook Programming

Leave a reply

The Blackberry Playbook tablet can be programmed, if you know the C language. I was quite surprised to discover how open it is. In case others are interested, here is the technique used:

– Download the Playbook software development kit (SDK) to your Windows PC. Obtain a code-signing key. That key can be used to publish apps to the app store. However, it is not necessary to publish an app in order to be able to use it on your own playbook.

– In the SDK, write your app. There are a number of sample apps available as code examples. For instance, I wanted to write a rotating torus application, to learn some graphics programming. Using the SDK, create a new project: File – New – Blackberry Tablet OS project, with the OpenGL ES 2.0 template. When built and run, that program creates a 2-D coloured square which rotates.

– Modify the template code. Strip out the square logic, and substitute code which defines the surface of a torus. The torus is formed from an N1-sided polygon, approximating a circle of radius R1, which is rotated in space, N2 copies, around a larger circle of radius R2. That forms the vertexes of the torus. Those vertexes are used to define triangular patches of surface, which in GLES2 are referred to as elements.

– Also include code which rotates the torus around the x, y, or z axes, modifying the rotating square code. Compile, debug, and run on the playbook until the code is working to satisfaction.

– BUT, a problem is encountered. The torus surface is not being drawn correctly. It turns out that the depth feature of GLES2 has not been enabled, in the screen initialization routine which is part of the bbutil.c library provided by Blackberry. It is necessary to modify bbutil.c to include a spec “EGL_DEPTH_SIZE, 24″ in the “attrib_list” specs within the bbutil_init_egl routine. Very fortunately, bbutil.c is provided as source code, and a distinct copy is set up with each newly created project, so it is easy to modify for this special situation. Now the depth testing works satisfactorily.

The result: an enjoyable coding project, and a nice rotating torus on the screen.

As a further refinement, code was added to detect the position of user touches on the screen and cause the torus to change its rotation, based upon the location of the screen touch. The test for the touch location is not as well documented as I needed, so here is the detail; Inside the “handleScreenEvent” routine, after a SCREEN_EVENT_MTOUCH… event has been identified, include a call to get SCREEN_PROPERTY_SOURCE_POSITION into a pair of short integers. Those will be the x and y coordinates of the screen touch, in pixels from the top left of screen.

Enjoy!

Ken Roberts

27-Nov-2013