Finding Munchausen Numbers with Perl

Boy has it been forever since I used Perl for something useful. It's been ages since I used Perl for anything actually. And the local Phoenix Perl Mongers group seems to have dissolved as well. That's too bad, because I really did enjoy going to those talks. But if nobody chips in for those talks the group quickly will die out. Eric Tank can only say so many things about the work he does and the things he knows. I wish I knew more about PDL to give a decent talk about it at a new Perl Mongers group meetings. 

But I digress. Today's blog post is a little blurb about Munchausen numbers. It is stated that in the base 10 number system there are only two known Munchausen numbers. I wrote a small Perl script to see if that was true for the numbers ranging form 1 to 500,000. I don't think I ever got to 500,000 before killing the script, but it was fun at least seeing one of the solutions that everyone already agrees upon. I learned some things about number theory. I would like to understand the proof behind the statement of Munchausen numbers. I currently don't understand it. Finding the Wikipedia article about the curious numbers killed my enthusiasm for the Perl script. Instead I had fun running the GAP script a paper provided into my old Debian Panasonic Toughbook So maybe later today I'll post the crappy Perl script I wrote along with the mysterious and amazing GAP code somebody else wrote. LInks and buttons to come in the near future folks! 

Playing with TrueBASIC

I was reading an old book about computer simulation of liquids and wrote this little script that uses the Lennard-Jones potential to calculate the random distribution of atoms in a gas lattice.

OPTION NOLET
! This code illustrates the calculation of the potential energy in a system of
! Lennard-Jones atoms.
! I create a set of atoms in random locations, calculate the Lennard-Jones potential,
! and then move the atoms in a way that would minimize this energy calculation.
N = 20    ! Number of particles
V = 0.0
EPSILON = 164
SIGSQ = .383
STEPX = 0.1
STEPY = STEPX
STEPZ = STEPY
STEP = 1
NOMOVE = 0
InitialPlotString$ = "IP = point(["
FinalPlotString$ = "FP = point(["
filename$ = "output.txt"
OPEN #1: NAME filename$, CREATE NEWOLD
ERASE #1


! Coordinate array for the particles
DIM RX(20)
DIM RY(20)
DIM RZ(20)

! Read coordinate data into arrays
PRINT "Initial atomic positions:"
PRINT "X","Y","Z"
PRINT #1: "Initial atomic positions:"
PRINT #1: "X","Y","Z"
FOR I = 1 TO N
    RX(I) = RND*10
    RY(I) = RND*10
    RZ(I) = RND*10
    PRINT RX(I),RY(I),RZ(I)
    PRINT #1:RX(I),RY(I),RZ(I)
    IF I <> N THEN
        InitialPlotString$ = InitialPlotString$ & "(" & STR$(RX(I)) & "," & STR$(RY(I)) & "," & STR$(RZ(I)) & "),"
    ELSE
        InitialPlotString$ = InitialPlotString$ & "(" & STR$(RX(I)) & "," & STR$(RY(I)) & "," & STR$(RZ(I)) & ")],color=red,size=20,opacity=1)"
    END IF
NEXT I

CALL LennardJones
PRINT "Initial V = "; V
PRINT #1: "Initial V = "; V
OLDV = V
DO WHILE (NOMOVE < 3*N)
    FOR I = 1 TO N
        ! Wiggle atom I in X and see if it reduces the potential.
        RX(I) = RX(I) + STEPX
!        PRINT "STEP ";STEP
        STEP = STEP + 1
        CALL LennardJones

        IF ABS(V) > ABS(OLDV) THEN
            RX(I) = RX(I) - 2*STEPX
!            PRINT "STEP ";STEP
            STEP = STEP + 1
            CALL LennardJones
            IF ABS(V) > ABS(OLDV) THEN
                NOMOVE = NOMOVE + 1
                RX(I) = RX(I) + STEPX
            ELSE
                OLDV = V
            END IF
        ELSE
            OLDV = V
        END IF

        ! Wiggle atom I in Y and see if it reduces the potential.
        RY(I) = RY(I) + STEPY
 !       PRINT "STEP ";STEP
        STEP = STEP + 1
        CALL LennardJones
        IF ABS(V) > ABS(OLDV) THEN
            RY(I) = RY(I) - 2*STEPY
!            PRINT "STEP ";STEP
            STEP = STEP + 1
            CALL LennardJones
            IF ABS(V) > ABS(OLDV) THEN
                NOMOVE = NOMOVE + 1
                RY(I) = RY(I) + STEPY
            ELSE
                OLDV = V
            END IF
        ELSE
            OLDV = V
        END IF

        ! Wiggle atom I in Z and see if it reduces the potential.
        RZ(I) = RZ(I) + STEPZ
        !PRINT "STEP ";STEP
        STEP = STEP + 1
        CALL LennardJones
        IF ABS(V) > ABS(OLDV) THEN
            RZ(I) = RZ(I) - 2*STEPZ
            !PRINT "STEP ";STEP
            STEP = STEP + 1
            CALL LennardJones
            IF ABS(V) > ABS(OLDV) THEN
                NOMOVE = NOMOVE + 1
                RZ(I) = RZ(I) + STEPZ
            ELSE
                OLDV = V
            END IF
        ELSE
            OLDV = V
        END IF
    NEXT I
    PRINT "STEP: ";STEP
LOOP

PRINT "Minimum V: ";OLDV
PRINT "Atomic Coordinates:"
PRINT "X","Y","Z"
PRINT #1: "Minimum V: ";OLDV
PRINT #1: "Atomic Coordinates:"
PRINT #1: "X","Y","Z"
FOR I = 1 TO N
    PRINT RX(I),RY(I),RZ(I)
    PRINT #1: RX(I),RY(I),RZ(I)
    IF I <> N THEN
        FinalPlotString$ = FinalPlotString$ & "(" & STR$(RX(I)) & "," & STR$(RY(I)) & "," & STR$(RZ(I)) & "),"
    ELSE
        FinalPlotString$ = FinalPlotString$ & "(" & STR$(RX(I)) & "," & STR$(RY(I)) & "," & STR$(RZ(I)) & ")], color=blue,size=20,opacity=1)"
    END IF
NEXT I
PRINT #1: "Atomic Positions for Sage plotting:"
PRINT #1: InitialPlotString$
PRINT #1: FinalPlotString$
PRINT #1: "Show(IP+FP)"
PRINT "Done."

SUB LennardJones
! This subroutine calculates the Lennard-Jones potential for the system
LOCAL I
V = 0.0
FOR I = 1 TO N - 1
    RXI = RX(I)
    RYI = RY(I)
    RZI = RZ(I)
    FOR J = I + 1 TO N
        RXIJ = RXI - RX(J)
        RYIJ = RYI - RY(J)
        RZIJ = RZI - RZ(J)

        RIJSQ = RXIJ^2 + RYIJ^2 + RZIJ^2
        SR2 = SIGSQ / RIJSQ
        SR6 = SR2 * SR2 * SR2
        SR12 = SR6^2
        V = V + SR12 - SR6
    NEXT J
NEXT I
V = 4.0 * EPSILON * V
PRINT #1: "STEP "; STEP; ":"
FOR I = 1 TO N
    PRINT #1: RX(I),RY(I),RZ(I)
NEXT I
!PRINT "V = "; V
PRINT #1:"V =";V
END SUB

END

Learning about CBM-BASIC

I've been solving Brilliant math problems with the simplest of BASIC programming languages lately. I do this because it's super convenient for me to write a quick little program to do some equation solving. CBM-BASIC is pretty nice for this because it works with Windows and Linux and my iPhone (by installing Hand-BASIC).

For a nice little example I solved this problem:

 With the help of Mathematica and CBM-BASIC. I used Mathematica to solve for z and plot the resulting equation to help me understand if there were a finite number of solutions. I was being lazy, because I could have figured this out just by plugging a few exploratory numbers in for X,Y, and Z to help me understand the bounds of the problem. Once I understood the bounds of the problem I was able to find all the solutions with CBM-BASIC in a relatively short time period.

Here is the resulting BASIC program:

10 REM Solve Brilliant Problem 2 of May 8, 2017 
20 S = 0
30 E1 = 0
40 MX = 150
50 MY = 50
60 MZ = 40
65 PRINT "FINDING SOLUTIONS..."
70 FOR X = 4 TO MX
80 FOR Y = 1 TO MY
90 FOR Z = 1 TO MZ
100 GOSUB 500
110 IF ABS(E1) < 0.1 THEN S=S+X+Y+Z : PRINT "S = ";S
120 NEXT Z : NEXT Y : NEXT X
130 PRINT "DONE." : END
500 REM EVALUATE FUNCTION E1
510 E1 = X^2*(Y^3+Z^3)-315*(X*Y*Z+7)
520 RETURN

And there you have it. I'll let the reader build CBM-BASIC and find the solution themselves. I'm having a lot of fun using this little Commordore 64 BASIC 2 implementation for my mental exercise. Who said BASIC causes brain damage?

Unicon benchmark results on Raspberry Pi Zero W

This was from a version of Unicon compiled from unicon-code-5082-trunk. A recent development snapshot.


Times reported below reflect averages over three executions.
Expect 2-20 minutes for suite to run to completion.

Word Size  Main Memory  C Compiler  clock    OS         
32 bit     424 MB       gcc 4.9.2   1.0 GHz  UNIX       

CPU                            
1x ARMv6-compatible processor rev 7 (v6l)                  

                                        Elapsed time h:m:s        |Concurrent |
benchmark                           |Sequential|   |Concurrent|   |Performance|
concord concord.dat                   01:0.685            N/A
deal 50000                            01:0.504            N/A
ipxref ipxref.dat                        35.286            N/A
queens 12                             01:13.523            N/A
rsg rsg.dat                           01:3.325            N/A
binary-trees 14                       01:36.725       04:55.517         0.327x
chameneos-redux 65000                      N/A       01:8.534
fannkuch 9                            01:8.319            N/A
fasta 250000                          01:22.580            N/A
k-nucleotide 150-thou.dat             03:19.441            N/A
mandelbrot 750                        08:10.612       13:54.547         0.587x
meteor-contest 600                    03:30.708            N/A
n-body 100000                         03:4.327            N/A
pidigits 7000                         14:37.462            N/A
regex-dna 700-thou.dat                01:52.136       03:33.752         0.524x
reverse-complement 15-mil.dat 
Run-time error 306
File reverse-complement.icn; Line 35
inadequate space in string region
Traceback:
   main(list_1 = [])
   getavgtimes(procedure run_reversecomplement,"15-mil.dat") from line 186 in run-benchmark.icn
   run_reversecomplement(list_5068944 = ["15-mil.dat"]) from line 129 in auxiliary.icn
   mapseq(">ONE Homo sapien...") from line 46 in reverse-complement.icn
   "
Run-time error 302
File reverse-complement.icn; Line 35
memory violation
Traceback:
   main(list_1 = [])
   getavgtimes(procedure run_reversecomplement,"15-mil.dat") from line 186 in run-benchmark.icn
   run_reversecomplement(list_5068944 = ["15-mil.dat"]) from line 129 in auxiliary.icn
   mapseq(">ONE Homo sapien...") from line 46 in reverse-complement.icn

Investigating Bayesian Network Modeling with Perl

I have a capstone project to finish this semester in order to get my Masters of Science in Engineering degree for Modeling and Simulation. It involves modeling a semiconductor fab using a combination of Bayesian network statistics and discrete event simulation. So far I'm in the literature review stage, although I do have a partially working discreet event simulation (DEVS) working from a previous course I took. So I was wondering how I'm going to run the Bayesian network analysis, and if there are any free open source packages already created for me to use. Well, it turns out that there is not much for Bayesian networks. Someone had asked this question on the Perl Monks forum Perl Monks forum but at that time they suggested the application was better suited using fuzzy logic or a naive-Bayes classifier. So per their advice I installed AI::FuzzyInference and AI::Categorizer::Learner::NaiveBayes onto my desktop. I have not had a chance to learn how to use them yet, but when I do I will be sure to document it here.

Check out Notebook Stories for a chance to win Monologue Notebooks!

My first application of the Perl Data Language

I've been installing and testing and playing with PDL for a few years now, but I've never had a need to use it for anything at work. Until now! Recently I needed to update the coordinate file for a set of devices to test. The layout editor showed that the centers of the devices were not rectangular shapes but polygons instead. That made it difficult for the software to give me a central coordinate of the object. What L-Edit does give me though are the vertices of the polygon that makes up that object, and that data was easy enough for me to copy and paste into TextPad. With a few regular expressions and search/replace commands I was able to format the data in a way that Perl could read. I then used PDL to quickly calculate the centroid of the coordinate set, and voila! I had my center coordinate for the device I needed to test. It was my first real use of PDL and I think it saved me a few lines of code. The cool thing was turning it into a PDL function, which could be called by the PDL REPL.

#!/usr/bin/env perl
use Modern::Perl;
use PDL;

sub FDC53_NE_9x9_X4Y0_Centroid() {
    my @V;
    push @V, [53028.301,22242.396];
    push @V, [53028.477,22241.163]; 
    push @V, [53028.532,22240.839];
    push @V, [53028.646,22240.268]; 
    push @V, [53028.958,22239.265];
    push @V, [53029.368,22238.308]; 
    push @V, [53029.672,22237.713];
    push @V, [53029.932,22237.280]; 
    push @V, [53030.687,22236.313];
    push @V, [53031.811,22235.189]; 
    push @V, [53032.783,22234.430];
    push @V, [53033.207,22234.176]; 
    push @V, [53033.814,22233.865];
    push @V, [53034.769,22233.456]; 
    push @V, [53035.763,22233.147];
    push @V, [53036.373,22233.025]; 
    push @V, [53036.628,22232.982];
    push @V, [53037.896,22232.801]; 
    push @V, [53038.608,22232.750];
    push @V, [53127.893,22232.750]; 
    push @V, [53128.603,22232.800];
    push @V, [53129.865,22232.981]; 
    push @V, [53130.134,22233.027];
    push @V, [53130.732,22233.147]; 
    push @V, [53131.735,22233.458];
    push @V, [53132.692,22233.868]; 
    push @V, [53133.287,22234.172];
    push @V, [53133.720,22234.432]; 
    push @V, [53134.687,22235.187];
    push @V, [53135.815,22236.315]; 
    push @V, [53136.566,22237.277];
    push @V, [53136.824,22237.707]; 
    push @V, [53137.134,22238.314];
    push @V, [53137.544,22239.270]; 
    push @V, [53137.852,22240.262];
    push @V, [53137.972,22240.861]; 
    push @V, [53138.020,22241.141];
    push @V, [53138.199,22242.396]; 
    push @V, [53138.250,22243.108];
    push @V, [53138.250,22332.893]; 
    push @V, [53138.200,22333.603];
    push @V, [53138.019,22334.865]; 
    push @V, [53137.973,22335.134];
    push @V, [53137.853,22335.732]; 
    push @V, [53137.542,22336.735];
    push @V, [53137.132,22337.692]; 
    push @V, [53136.828,22338.287];
    push @V, [53136.568,22338.720]; 
    push @V, [53135.813,22339.687];
    push @V, [53134.681,22340.819]; 
    push @V, [53133.727,22341.564];
    push @V, [53133.293,22341.824]; 
    push @V, [53132.686,22342.135];
    push @V, [53131.731,22342.544]; 
    push @V, [53130.737,22342.853];
    push @V, [53130.127,22342.975]; 
    push @V, [53129.872,22343.018];
    push @V, [53128.604,22343.199]; 
    push @V, [53127.892,22343.250];
    push @V, [53038.607,22343.250]; 
    push @V, [53037.897,22343.200];
    push @V, [53036.635,22343.019]; 
    push @V, [53036.366,22342.973];
    push @V, [53035.768,22342.853]; 
    push @V, [53034.765,22342.542];
    push @V, [53033.808,22342.132]; 
    push @V, [53033.213,22341.828];
    push @V, [53032.780,22341.568]; 
    push @V, [53031.813,22340.813];
    push @V, [53030.686,22339.686]; 
    push @V, [53029.933,22338.722];
    push @V, [53029.676,22338.293]; 
    push @V, [53029.365,22337.686];
    push @V, [53028.956,22336.731]; 
    push @V, [53028.647,22335.737];
    push @V, [53028.525,22335.127]; 
    push @V, [53028.482,22334.872];
    push @V, [53028.301,22333.604]; 
    push @V, [53028.250,22332.892];
    push @V, [53028.250,22243.108]; 
   
    my $data = pdl(@V);
    my $mySums = sumover $data->xchg(0,1);
    my @dims = $data->dims;
    return $mySums / $dims[1];
}
       
1;

Converting Excel Macros to Perl

I'm currently working on porting our Excel macros at work to LabView. To get a quick prototype up and running I have been translating the Visual Basic code and the formulas in all of the spreadsheet cell to Perl 6. In the process I get to learn more about Perl 6 and see what all the fuss is about. The Excel macro takes a directory of data files from our transistor testing procedure and extracts certain model parameters for our design engineers. The data is located in tab separated plain text files. They contain voltage and current data for different aspects of the transistor. Each folder contains the data of over a dozen transistors. Here is what I've done so far:

#!/usr/bin/env perl6
use v6;

my $Append = False;
my $Test_Descriptor_Keys_File = "(compute)";
my $Lot_ID = "E1517-002";
my @Wafer_List = "Wafer_11", "Wafer_12";
my @Analysis_Master_Name = "AnalysisMaster_Autoprobe_Prod_20111005";
my $Test_Name_For_Processing = "COW";
my $Path_To_Data_File = "./";
my $Process_Reverse_Sweeps = 1;  my $Print_Graphs = 1;
my $Print_Graphs_To_File = 1;
my $Delete_Summary_Filename_Sheets = 0;
my %Environment = Lot_ID => $Lot_ID, Wafer_List => @Wafer_List, Path_To_Data_File => $Path_To_Data_File,
    Test_Name_For_Processing => $Test_Name_For_Processing, Analysis_Master_Name => $Analysis_Master_Name,
    Test_Descriptor_Keys_File => $Test_Descriptor_Keys_File, Append => $Append, Process_Reverse_Sweeps => $Process_Reverse_Sweeps,
    Print_Graphs => $Print_Graphs, Print_Graphs_To_File => $Print_Graphs_To_File, Delete_Summary_Filenames_Sheets => $Delete_Summary_Filename_Sheets;

sub main() {
    say "Starting...";
    my $myReturn = 1;
    my $myInput = "";
    my $PrintGraphsStatus = "(yes)";
    my $PrintGraphsToFileStatus = "(yes)";
    my $ProcessReverseSweeps = "(yes)";
    my $AppendDataToSummary = "(no)";

    while ($myReturn) {
        say "Choose Function: ";
        say "1. Process TFT Autoprobe Data";
        say "2. Process TFT Princeton Data";
        say "3. Print Graphs? " ~ $PrintGraphsStatus;
        say "4. Print Graphs to File? " ~ $PrintGraphsToFileStatus;
        say "5. Process TFT Data";
        say "6. Delete Summary & Filenames Sheets";
        say "7. Process Reverse Sweeps " ~ $ProcessReverseSweeps;
        say "8. Append Data to Summary " ~ $AppendDataToSummary;
        say "9. Quit";

        $myInput = get();

        given ($myInput) {
            when 1 { $myReturn = &Process_TFT_Autoprobe_Data(%Environment); }
            when 2 { $myReturn = &Process_TFT_Princeton_Data(%Environment); }
            when 3 { $myReturn = &Print_Graphs(%Environment); }
            when 4 { $myReturn = &Print_Graphs_To_File(%Environment); }
            when 5 { $myReturn = &Process_TFT_Data(%Environment); }
            when 6 { $myReturn = &Delete_Summary_Filenames_Sheets(%Environment); }
            when 7 { $myReturn = &Process_Reverse_Sweeps(%Environment); }
            when 8 { $myReturn = &Append_Data_to_Summary(%Environment); }
            when 9 { exit(0); }
            default { $myReturn = 1; }
        }
    }
}

sub Process_TFT_Autoprobe_Data() {
    my %Environment = shift;
    my $Lot_ID = %Environment{'Lot_ID'};
    my @Wafer_List = %Environment{'Wafer_List'};
    my $Path_To_Data_File = %Environement{'Path_To_Data_File'};

    say "Enter your lot ID:";
    $Lot_ID = get();



    say "called Process TFT Autoprobe Data";
}

sub Process_TFT_Princeton_Data() {
    my %Environment = shift;
    say "called Process TFT Princeton Data";
}

sub Print_Graphs() {
    my %Environment = shift;
    say "called Print Graphs";
}

sub Print_Graphs_To_File() {
    my %Environment = shift;
    say "called Print Graphs to File";
}

sub Process_TFT_Data() {
    my %Environment = shift;
    say "called Process TFT Data";
}

sub Delete_Summary_Filenames_Sheets() {
    my %Environment = shift;
    say "called Delete Summary Filenames Sheets";
}

sub Process_Reverse_Sweeps() {
    my %Environment = shift;
    say "called Process Reverse Sweeps";
}

sub Append_Data_to_Summary() {
    my %Environment = shift;
    say "called Append Data to Summary";
}

main();
1;

Whole Cell Simulation Project

Lately I've been thinking about computer simulations and what would be the most fruitful thing to simulate right now. I've read about simulating an entire bacterial cell at the molecular lever in the past, but only recently have I learned that the people who accomplished that feat have made their software available for free downloading. I was super excited and quickly installed the source code into three computers. The source code is found here: https://github.com/CovertLab/WholeCell

The code is written in Matlab 2010b. I've been successful at getting it to work on Matlab 2015a and Matlab 2010a, with Windows 7 and Redhat Linux respectively. Oddly I have not been able to get the software to work on Scientific Linux 6 running either Matlab 2012a or 2010b. That machine runs a dual core Pentium D processor and there seems to be something wrong at a deep level with one of the number crunching libraries.

Quick and dirty Matlab/Octave replacement using JScript

In my previous post I showed that it is possible to make a simple read-eval-print-loop (REPL) in JScript for the Windows command line. Today I show that it is also possible to add libraries to this REPL for quick access and experimentation. I have been adding functions to my JScript matrix library for work when I realized that I was quickly needing Matlab levels of functionality in order to get what I needed. Today I decided to try loading a Javascript library for numerical work into JScript. I was surprised how painless the process was!

The first thing I needed to do was write a WSF file to load the library with the jshell REPL I created in my last post. I then found a nice Javascript library online (http://numericjs.com/numeric/index.php) and integrated this library with my jshell.js REPL script into the WSF file. I called my file jlab.wsf.