vsoch

Madlib Python Game!

In my free time I’ve been figuring out object oriented programming. I am happy to post that I have just finished my first OOC (object oriented code) in python, and am excited to share what I’ve learned!

I decided to pick a simple project that would require some file reading and writing, as well as user input and standard output, since these basic ideas can be tweaked and applied towards running imaging analysis. The following script uses a text file as a template to play a game of madlib. The text file can be however you like, however the words that will be found and presented to the user for substitution should be in the format SUB_NOUN_SUB, SUB_ADJECTIVE_SUB, etc.

Here is how it works:

1. When madlib.py is called, we of course start in main(), and the madlib class makes an instance of itself (mymadlib = madlib())
2. mymadlib.getFile(): The user is asked to specify a text file to be used. If in the same directory, you can just type “mytext.txt.” If not in the same directory, the full path should be written out. The script checks to make sure that the file is a .txt or .doc file (other types could be added), and then checks if the file exists. If it does, we open it, if not, we spit out an error.
3. mymadlib.outputName(): We then keep asking the user for an output name until one is entered that doesn’t already exist in the current directory. The output will be put in the present working directory.
4. mymadlib.readFile(): This function does the majority of the heavy lifting. It first reads all words in the text file into a list, and and then uses regular expressions to create an expression to match of the format (SUB_anything_SUB). Case does not matter, however there can be no spaces between the SUBs. If there is a space, the expression will get split into two when the file is read in, and then will not get matched, and will be not filled in with user input. So it’s important to format the substitutions with no spaces! Underscores are fine.
5. We then check if the length of the list is greater than zero, meaning we aren’t dealing with a blank file. When I tested this, it also reported a file filled with only carriage returns as having a length of 0,m so it should skip over reading the words if the file is empty OR filled with carriage returns.
6. We search through each word looking for the expression. If we find a word to substitute, we present it to the user with
7. **mymadlib.getWord(): **which is referred to as self.getWord() since we are calling it within the class. This function simply prompts the user for input for whatever expression is specified between the SUB’s.
8. It returns the user’s word choice, which is then insert into the list at the same spot as the original word, and then since the original word is moved one spot forward in the list, we remove it completely by using a function called pop() which can be applied to lists.
9. When this process has been done for every SUB word, we use **mymadlib.createLib() **to print out new list to the terminal for the user to say, and also save it to our specified output file. Since this function is also called within the class, we reference it as self.createLib().

A few things that I learned:

• A class is essentially an object. You create instances of the object, and then can call functions specific to that object on your instance. So for example, since my class is called class madlib:, I can create a madlib object with a simple mymadlibname = madlib(). Then if I want to call a function in this class on my object instance, I can say mymadlibname.myfunction().
• If we create an object and simply call the class as a script (calling madlib.py, for example) this goes directly to the main() class.
• The if statement at the bottom is necessary because it makes sure that main() is only run if the script is being called. If madlib is just being used as an object for another script, then we don’t want to call the first main().
• The import statements at the top are basically pointing to common classes that add functionality to my script. For example, import re let’s me use regular expressions. How do I feel about the ability to do this and select from a seemingly endless library of classes? …like a kid in a candy store! :O)
• Each function within the class, since we are dealing with an object, by default takes the object as an input parameter, represented as ‘self.’ If a function has multiple inputs ‘self’ MUST be the first parameter, unless it’s called with the variables specified, in which case the order doesn’t matter.
• A couple of very powerful data structures include lists, dictionaries, and tuples. I created a tuple by accident, and had a lot of fun playing around with these different structures with my friend!
• re.compile returns a MatchObject, which can be referenced with myobject.group() (to print out the entire match) and myobject.start() and myobject.end() to print the start and end locations of the match in the string.
• Holy cow, if I just type “python” into my terminal window, that allows me to play! Previously I was writing python code in a .py file, and testing it by running it.

Here is the class and a sample text file, if you would like to try out my script! If you make any madlib’s, send them my way and I’d love to play them :O)

Madlib.py
TownStory.txt

Beer Rankings 1.0 (beerrank.m)

I should say first and foremost that I do not drink! However, my lab is going into the second year of what will be an annual Christmas tradition of having a “Holiday Beer Tasting,” and I am the official server. Last year we recorded everything on paper and had to do all the calculations manually, but this year I decided it would be cool to go a little more high tech. I of course would never pass up the opportunity to write a script!

Since yesterday evening I’ve whipped up and tested beerrank.m, which will take in all beer rankings, calculate averages, record raw and average data to file, and present the winners. The script works as follows:

**1) Run **by typing beerrank into the matlab terminal window.

**2) Input **the number of beers, and the desired name of the output file. If the output file already exists, the data will be appended to the next empty row. If the file does not exist, it will be created with all the appropriate headers (what I recommend to do!) The file will be put in whatever is your pwd (present working directory).

3) **Input Beer Names:** A window will pop up for each of the number of beers that you have specified! In this window you should put the actual beer name. Note that this step should be done by the server whom knows the secret identify of each beer in the order that it is being served. These names will be hidden from view until the final results are displayed, and beers will be referenced as BEER 1, BEER 2, etc.

**4) **Next, the user is prompted with an empty rankings table for BEER1. The ranking categories have been determined in advance, and at this point the server can start serving, and one person can input results as they are decided by each person, or you could go around the table at the end and have everyone read the ranking for each category. Note: An update to the script would be to allow the user to specify categories.

**5) Here **is an example of input data. Numbers can be entered with spaces between them, and extra spaces on the front, end, and in between do not matter – the script will nip them in the butt! Note that there is currently no input validation – so in the case that a user entered a non numerical value anywhere other than the comments box, the script would error out when it converts the string number into a double (a type of number) with the function str2double(). This could be another future update – to save data temporarily, check it, and re-present the prompt to the user in the case of a faulty input. I also want to note that it doesn’t matter how many numbers are input into each ranking – the average will be calculated nonetheless. I thought it would be unwise to ask the user to specify a number of rankings in advance, and then only accept that number!

6) Output files **produced include a **RESULTS.mat, a matrix of all the variables, which is updated as the script runs in the case that someone accidentally shuts out of it, and two excel files: one for raw data, and one for the averages.

**7) **After the output files are produced, the script calculates the winners for each category, and presents them to the user, along with the revealed list of beer names. Note that I didn’t do any sort of rounding or truncation for decimal points.

**8 ) Average Results File **contains the average ratings for each category for each beer, in the case that people are curious about 2nd, 3rd, 4th place, etc. Here is also where the comments are stored, which could be fun to save to look back on every year, if you’re into that kind of thing. I should also point out that the script doesn’t do anything in the case of two beers coming in at a tie – it simply keeps whichever one was found first as a winner. This is another update that could be implemented at some future point.

**9) Raw Data Results File: **I’m a proponent of always having your raw data somewhere, in the case that there was an error in calculation, or the results file explodes, you still have something to work with! Here it is.

and here is the script, if you care to take a look at the hard-coded monstrosity, haha. It’s shamefully ugly, but it’s going to work perfectly for how we need it, and I’m sure I’ll make a second, more elegant version at some future point!

Merry Christmas everyone! Be safe!

**Added December 15, 2010: **Here are the actual results:

and for lab pictures, see here

Simple Fudge

During a busy Christmas-time, I find myself in an incentive conundrum. On the one hand, there is a lot of work to do, and projects that I want to be working on when the work dwindles down. On the other hand, I want to devote some time to making something nice for Christmas gifts. What is the solution to this predicament? Fudge, of course! It has the highest deliciousness / time-spent-making ratio, and is a relatively affordable investment since a small number of ingredients go very far. There are also no constraints on creativity, meaning no limit to what you can dream up, combination wise. Lastly, it travels and packages well, and it’s a quick, no-bake delicacy, meaning that you can whip up many batches with only a pot and a pan or two. For these reasons, I’ve made different variations of fudge for years, and it never gets old!

Tools Required:

• a medium sized pot
• stove
• a mixing spoon
• a square pan (8X8 or 9X9 works well) – square edges work better
• aluminum foil (optional)
• a spatula (optional)

Ingredients Required:

• 3 cups chocolate bits, which is 1 1/2 bags (I prefer semi-sweet)
• 1 can sweetened condensed milk
• a smatter of salt
• about 1 – 1/2 tsp vanilla extract

Instructions:

**1. **Find a friend to cook with. Mr. Kitchen Timer will have to do!

**2. **Obtain ingredients: chocolate chips, vanilla, salt, sweetened condensed milk

**3. **Obtain your mix in of choice (think candy, nuts, cookies, marshmallows, etc) and prepare it appropriately. For large items like nuts and cookies, you likely want to cut the mix-in into smaller pieces. For oreos I like to cut them into quarters. Finish cutting your oreos, and set aside 16-25 quarters for top-decoration.

**4. **Set the rest of the oreos aside, and prepare a pan by throwing a piece of aluminum foil in it. Be sure to have edges peeking out so when the fudge is set, you can use them to easily lift the block out of the pan.

**5. **Next, add your 1 + 1/2 bags chocolate bits, 1 can sweetened condensed milk, and smatter of salt to a medium pot, and set on the stove on low. This is also where you should assess the ratio of chocolate to sweetened condensed milk, depending on your mix-in. Subtracting condensed milk or adding chocolate bits will make the fudge more sturdy and dense, which might be good for runnier or floopy mix-ins, but you lose some of the softness and shine. For 0reos, I like to make a little denser fudge to match the cookie, but for most other things I stick to the normal ratio. Don’t forget te salt – it is one of the essential ingredients to making the chocolate flavor really come out!

**6. **You are going to want to start mixing almost immediately. The ingredients will melt together very quickly, and a little arm action will evenly mix everything and prevent the bottom from burning. Keep mixing as you watch the fudge go through four stages: separate and cold, coated, melting and lumpy, and finally, uniform.

**7. **Just when the fudge is uniform (and not any later!) you want to remove it from the heat, and mix in the vanilla quickly. If you add the vanilla before the end point I’m pretty sure that it will cook off! If you are making a fudge with a normal ratio of chocolate / sweetened condensed milk, it will look shiny (see batch below, compared to batch above).

**8. **This is where you need to work quickly. Add the mix-in, mix it just enough to evenly distribute, but not so much that the mix-in starts to break down. Then, plop into your prepared pan, and if you have a denser fudge, use plastic baggies to press it down flat. If you have a creamy, shiny fudge, it should be easily spreadable with a spoon or spatula.

**9. **I would only recommend the “squishing” technique shown above if you plan to cover it with something else. Bare, it’s a little fugly. You next want to add your toppings, being sure to press in slightly so they stick, and then put in the refrigerator to chill. Do NOT put in the freezer – you will regret it! You can see below that I had already made a peppermint white chocolate version. I should note that white chocolate is a very fussy fudge base to work with. If you add a normal amount of sweetened condensed milk it will have the consistency of peanut butter, and never harden. The trick is to use more white chips, and add a little margarine and just enough sweetened condensed milk so it mixes and melts. You will want to get it off the stove as soon as it melts, otherwise it burns almost immediately. Then you have to work extremely quickly to get it into the pan, and you definitely will be pressing as opposed to spreading. But the deliciousness of white chocolate and peppermint is worth it!

**10. **The nice thing about fudge is that you can work on other batches while your previous ones are setting. It takes usually an hour to two, depending on the consistency. For my last batch, I decided to combine white and regular chocolate chunks, and use Almond Joy pieces as my mix-in.

**11. Note **how the fudge is shiny and easily spreads into the pan, with a pretty texture on the top. This is all about the ratio of chocolate < > sweetened condensed milk.

**12. **Place with the other pans to chill. Note that I have used the hanging aluminum foil to lift the center block from the pan, since I needed to use it for the new batch. This means that it is time to cut into squares! First, you want to hold the block in one hand, and peel away the wrapper from the bottom (as much as you can with your hand there). Then, you want to use a long, sharp knife to cut the fudge. I like to cut off the edges to give each piece a nice, clean cross section, but you don’t have to. Your cutting strategy will vary based on the consistency and mix ins. I generally find that pressing and gentle rocking, and then using force to “chop” it works best. Smaller squares are usually preferable to consumers, and it always makes your batch go a little farther. Lastly, put them into paper wrappers or a container.

**13. **This is a good time to wash dishes. If you didn’t do various colorings, you should mainly just have the pot, a spoon, can-opener, and your knife.

**14. **For packaging, I like to make a variety, and put in cute little bags. You can do one or two days in advance, as long as you keep them refrigerated, and sealed air-tight. These are for my lab-mates!

15. Merry Christmas!

My Gum Chewing Algorithm :P

I was writing to a friend and somehow got into the topic of chewing gum, something that I’m very fond of. Matter of fact, I have a very methodical way of consuming my favorite three flavors, Orbit Cinnamon, Trident Strawberry, and Extra Watermelon, while I’m working busily away. So instead of writing it out in words, I decided to write a script in faux code. It was fun! And mostly for it’s preservation so that I can find it in X years down the road, here it is!

% We start out with zero pieces of gum, and nothing being chewed.

% The order keeps track of my special order of chewing, to be cycled through, 1 through 5
pieceCount = 0;
piecesLeft = 7; % we start with 7 pieces
chewing={}; % we start chewing nothing
order=1; % start with order 1
while (piecesLeft > 0)
% Always spit out gum when we have 2 pieces that have lost flavor:
if ((pieceCount == 2) && (lostflavor(chewing)==’true’)
spit(chewing);
% Or spit if we have one extra watermelon piece that has lost flavor
elseif((pieceCount == 1) && (lostflavor(chewing)==’true) && (chewing==’extra(1)’||’extra(2)’))
spit(chewing);
else
order = consume(order);
end
% Here is the function to chew the gum based on our current order
function consume(order_number)
switch order:
% Start with one strawberry piece, and chew until it’s no good
case 1: eat(strawberry(1))
while(lostflavor(strawberry(1))~=’true’)
chew()
end
% When it loses flavor, add an extra watermelon piece
add(extra(1),strawberry(1))
while(lostflavor(strawberry(1),extra(1))~=’true’)
chew()
end
return 2;
% Then eat the last extra piece, and chew until it’s no good
case 2: eat(extra(2))
while(lostflavor(extra(2))~=’true’)
chew()
end
return 3;
% Eat both strawberry pieces together, because they are small!
case 3: eat({strawberry(2),strawberry(3)})
while(lostflavor({strawberry(2),strawberry(3)})~=’true’)
chew()
end
return 4;
% Grand finale is two cinnamon pieces!
case 4: eat({orbit(1),orbit(2))
while(lostflavor({orbit(1),orbit(2)})~=’true’)
chew()
end
return 5;
case 5;
fprintf(‘%s\n’,’We’ve finished all of our gum! Replenish supply and start over’)
exit;
end
end
% Function to eat gum
function eat({gum})
nomnom(gum);
piecesLeft=piecesLeft-size(gum);
chewing=gum;
end
% Function to spit out gum being chewed
function spit(tospit)
trash(tospit)
pieceCount = 0;
piecesLeft = piecesLeft -size(tospit);
chewing = {};
end
end
It kind of makes me want to try writing scripts for little things in life that are methodical like brushing your teeth, cooking, or running. It also makes me realize that people who design the machines and gadgets that we use in our everyday life actually DO think about these sorts of algorithms. That’s so cool! :O)

Loading Dynamic Variables in a Static Workspace in MATLAB

I’m working on a set of scripts to do our full resting connectivity analysis with the conn box from MIT, and of course this involves creating and processing design matrices on the cluster, and then when all is said and done, going back and changing the cluster paths to local ones. The scripts themselves are in a testing phase (and I will probably write about them when they are more finalized) but I wanted to share a silly error that I encountered while working on my script that changes paths in the design matrix.

The connectivity box produces an output .mat file (the design matrix that can be opened again in the GUI) that has a user specified name, so obviously the name will be stored within a variable in the script. Step 1 in my path changing script is to load this matrix. In the MATLAB GUI this is a simple load(mat_name), where mat_name might be something like “conn_rest.mat.” However, when I tried to do this in my script, I got the following error:

??? Error using ==> load
Attempt to add “CONN_x” to a static workspace.
See MATLAB Programming, Restrictions on Assigning to Variables for details.
Error in ==> conn_change_paths>change_conn at 151
load(mat_name’)
Error in ==> conn_change_paths at 136
change_conn(mat_name,oldpath,newpath,slash,old_slash);

I would guess that this is a common error for many MATLAB users granted that many might want to load .mat files or variables that may have dynamic names in a static context at one time or another. What I don’t quite understand is what is so dangerous about doing this? The variable that I wanted to load from the .mat, called CONN_x, will always have the same name and fields. Perhaps MATLAB has a fear of loading mysterious matrices with an unknown number of variables? So, for my next attempt, I tried to reassure MATLAB that I only wanted to load one variable, and it had a static name (CONN_x):

load(mat_name,’CONN_x’)

which is telling it specifically to load the CONN_x variable from the .mat that is specified. I got the same error. So it seems that the only solution is to load CONN_x into a structure, doing something like this:

C= load(mat_name,’CONN_x’)

and of course this means that my variable is now references as C.CONN_x.fieldname instead of just CONN_x.fieldname. Given that the connectivity toolbox is hard-coded expecting CONN_x, we can’t have this extra C hanging around. We can fix this problem with:

CONN_x = C.CONN_x

The entire thing seems sort of silly, because MATLAB allows me to do something likeload(‘REX.mat’,’params’)in a static context but not the exact same command using a variable .mat name. Matter of fact, and I can’t believe that I’m writing about something this silly! I decided to because I can imagine that many users would encounter this error, and in the rare case that a frustrated person stumbles on this post, it is well worth it!

Journal Club Discussion Functional Connectivity in ASDs

I gave my first presentation in about three years to present a paper for Journal Club. It by no means was perfect (and I apologize to Dr. Cohen for calling him “Sasha” instead of “Simon”), but it was a lot of fun! What I learned is that there is a lot of value in just presenting or talking about a paper with a group of people, as each person brings an individual knowledge set, questions, and opinions. It can lead to ideas for new research, and a better understanding of the discussed paper than is obtained from reading the paper in isolation.

Part I

Part II

Part III

End

Grid Matrix!

I’ve been busy with projects for the lab, but I got the chance to do one quick coding challenge this weekend. The challenge presented a 20 X 20 matrix of numbers, and asked:

What is the greatest product of four adjacent numbers in any direction (up, down, left, right, or diagonally) in the 20×20 grid?

As an example, the numbers 26, 63, 78, and 14 were highlighted in red.

I particularly like this challenge because it’s very cut and dry. I know that if I create an algorithm that goes through all the possibilities, I will undoubtedly find the answer! And I had to use matlab for this one… because it is the Matrix Laboratory! It was a nice little bit of fun after a Sunday filled with running analysis and practicing a journal club presentation. Ok, I know, I shouldn’t try to deceive anyone that I suffered though that. I very much enjoyed both those things, otherwise I wouldn’t have chosen to do them!

function grid_product()
%--------------------------------------------------------------------------
% This function finds the greatest product for any four diagonal, vertical,
% or horizontal numbers in a 20 X 20 grid
%--------------------------------------------------------------------------

% Read in the grid matrix
grid_matrix ={08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08; 49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00; 81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65; 52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91; 22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80; 24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50; 32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70; 67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21; 24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72; 21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95; 78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92; 16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57; 86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58; 19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40; 04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66; 88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69; 04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36; 20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16; 20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54; 01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48;};

highestproduct = 0;

% DIAGONAL RIGHT: Go through the grid, and for each spot, find the product of the diagonal.
for i=1:17

for j=1:17
product = grid_matrix{i,j}*grid_matrix{i+1,j+1}*grid_matrix{i+2,j+2}*grid_matrix{i+3,j+3};
if (product > highestproduct)
highestproduct = product;
end
end
end

% RIGHT and LEFT: Go through the grid, and for each spot, find the % product of horizontal numbers.
for i=1:20
for j=1:17
product = grid_matrix{i,j}*grid_matrix{i,j+1}*grid_matrix{i,j+2}*grid_matrix{i,j+3};
if (product > highestproduct)
highestproduct = product;
end
end
end

% DIAGONAL LEFT: Go through the grid, and for each spot, find the product of the diagonal.
for i=1:17
for j=4:20
product = grid_matrix{i,j}*grid_matrix{i+1,j-1}*grid_matrix{i+2,j-2}*grid_matrix{i+3,j-3};
if (product > highestproduct)
highestproduct = product;
end
end
end

% UP AND DOWN: Go through the grid, and for each spot, find the vertical
% product,vstarting at the top.
for i=1:17
for j=1:20
% Check if the farthest out value exists, 4 rows down, 1 to the
% right
product = grid_matrix{i,j}*grid_matrix{i+1,j}*grid_matrix{i+2,j}*grid_matrix{i+3,j};
if (product > highestproduct)
highestproduct = product;
end
end
end

fprintf('%s%d\n','The highest product is ',highestproduct);
end

More Number Games

These are incredibly fun, and I realize that it doesn’t make sense to write a separate entry for each one, so I will post them in clusters. I decided to attempt most of these using java over matlab, which will be good review, and more challenging than Matlab. Keep in mind that I haven’t looked at java in 6 years, so the first few will likely be done by brute force! For the ones done in Matlab, this is because I was working on a computer other than my personal laptop.

Smallest Number Divisible by 1 though 20
Pythagorean Triplet (Matlab)
Count Primes (Matlab)
Fibonacci Even Sum
Difference b/w Sum of Square and Square of Sum
Largest Palindrome Produced by Multiplying Two 3 Digit Numbers
Sum Multiples of 3 and 5 Below a Specified Ceiling (Matlab)

Smallest Number Divisible by 1 through 20

This challenge was to find the smallest positive number that is evenly divisible by all of the numbers from 1 to 20. Surprise! I didn’t use a script. I found this question refreshing because I was able to whip out pencil and paper, and solve it that way. I first listed the numbers from 1 and 20, and the prime factors of each:

So, it makes sense that, if you have a number like 16 which breaks down to 2222, this also includes the prime factors of 2 (2) and 8 (222), so since we are finding the smallest number that is divisible by the numbers 1 through 20, I would only need to include 16 (2222) in this calculation, and that will by default be divisible by 2 and 8 as well. So, when I had these listed out on paper, I just circled the largest occurrence of each prime factor, and mutiplied these together to get the answer. I won’t put the answer here, but the numbers I multiplied were:

(2222)(33)(5)(7)(11)(13)(17)*(19)

Pythagorean Triplet

function pythagorean_triplet(number)
%————————————————————————–
% This function finds a pythagorean triplet that adds up to a certain
% number input by the user (number), if one exists!
%
% I solved the two equations, both for (c*c), and then set them equal
% to one another, and cycled through potential sets of a's and b's
% until a solution was found that sets both sides of the equation equal to
% one another, the value of (c*c)
%
% EQUATION 1:  a + b + c = 1000    –> so EQUATION 3: c = 1000-a-b
% EQUATION 2: ((a*a)+(b*b))=(c*c)  –>; Pythagorean Theorum
%
% Square EQUATION 3    (c*c) = ((1000-a-b))*(1000-a-b))
%
% Set EQUATION 2 and EQUATION 3 equal to one another, since both
% equal (c*c).  Numbers a and b that satisfy that equation = solution!
% ((a*a)+(b*b)) = ((1000*(-a)*(-b))(1000*(-a)*(-b)))
%————————————————————————-

for a=1:number
for b=1:number
if (((a*a)+(b*b)) == ((number -a -b))*(number -a -b))
fprintf('%s%d\n','a = ',a);
fprintf('%s%d\n','b = ',b);
fprintf('%s%d\n','c = ',(1000-a-b));
fprintf('%s%d\n','The product is: ',(a*b*(number-a-b)));
return;
end
end
end

fprintf('%s%d\n','There is no pythagorean triplet that adds up to ',number);
end

</code> </pre>

Count Primes

This script returns the nth prime number, where n is specified by the user. Since I had already created a function that checks if a number is prime for a previous challenge, I decided to modify that script for this task. The script sets a counting value, j, equal to 1, and then enters a while loop that continues until the variable ‘done’ is set to yes. Within this while loop for every odd number, we check if the number is prime. If the number is prime, we add one to the count of the prime_counter variable, which keeps track of what prime number we are currently at. When this prime_counter is equal to the user specified number, this is our solution.

function prime_count(number)
%--------------------------------------------------------------------------
% This function finds the nth prime number, where n = a number input
% by the user
%--------------------------------------------------------------------------

% Place the original number in a variable so we can reference it later
prime_counter = 0;

% First check if we have been given 0 or 1.
if (number == 0)
fprintf('%d%s\n',number,' is not a valid input number.');
return
end

if (number < 0)
fprintf('Please enter a positive number.\n');
return
end

done = 'no';
j = 1;

while strcmp(done,'no')
%If the number is even, skip it
if (mod(j,2)==0)
else

% Figure out if we have a prime
if (isprime(j))
prime_counter = prime_counter+1;
if prime_counter == number;
fprintf('%s%d%s%d%s\n','The prime number at the ',number,' spot is ',j,' and we have finished searching for primes!')
done = 'yes';
end
end
end
j = j+1;
end

%--------------------------------------------------------------------------
% Function isprime
%
% This function checks if an input number is prime by dividing by all
% possible factors through the input number / 2.  In the case that there
% is always a remainder for all these divisions, we know that the number
% has to be prime, meaning that it is only divisible by 1 and itself.
%--------------------------------------------------------------------------
function f = isprime(number_two)

% if our number is even, we don't continue
if mod(number_two,2)==0
f = 0;
return;
else
for i = 3:(number_two/2)
if mod(number_two,i) == 0
f = 0;
return;
end
end
f = 1;
return
end
end
end

Fibonacci Even Sum

This was the first challenge that I attempted with java. I literally downloaded netbeans right before writing this, and reintroduced myself to java syntax, which I haven’t touched in 6 years! It was super awesome, and came back to me pretty quickly. While I’m still very rusty, I apologize for the lack of elegance… I will need time to warm up!

This script calculates the numbers in the fibonacci sequence up to a certain ceiling, specified by the user. It first checks to make sure that there is only one input argument, and that it is an integer. In then calculates the Fibonnaci numbers up to that number, and in the case that the number is even, adds it to a sum. When we reach the ceiling, the script returns the sum to the user. It was very cool to compile and build this, so that I could run it in the command prompt!

/*
* Fibonacci Even Sum
* This script finds the sum of all even Fibonacci numbers below a user specified number (4 million)
*/
package fibevensum;

/**
* Vanessa Sochat
* October 23, 2010
*/
public class Main {

public static void main(String[] args) {

// Check to make sure we have one arguments
if (args.length != 1) {
System.out.println("This script only takes one input argument.");
System.out.println("...a number ceiling to check for primes below it.");
System.exit(0);
}

//Check to make sure we have an integer
try {
int ceiling = Integer.parseInt(args[0]);
System.out.println("The number we are using is " + ceiling);
} catch(NumberFormatException ife) {
System.out.println(args[1] + "is not a number, exiting");
System.exit(0);
}

// Start off by declaring the first two numbers, and a sum.
int one = 1;
int two = 1;
int holderone;
int holdertwo;
int sum = 0;

// While our lower number (one) is below the ceiling, keep going
while (one < Integer.parseInt(args[0]))
{

// If both numbers are below the ceiling
if (one < Integer.parseInt(args[0]) && (two < Integer.parseInt(args[0])))
{

if (one==0)
{
System.out.println("Adding " + one + " to the sum");
sum = sum + one;
}

if (two%2==0)
{
System.out.println("Adding " + two + " to the sum");
sum = sum + two;
}
}

// If the smaller number is below the ceiling, but
// the larger is not
else if (two >= Integer.parseInt(args[0]) && one < Integer.parseInt(args[0]))
{
if (one==0)
{
System.out.println("Adding " + one + " to the sum");
sum = sum + one;
System.out.println ("The sum is " + sum);
System.exit(0);
}
}

// If BOTH numbers are above the ceiling (although we shouldn't be
// in this loop at all) return the sum)
else if (two >= Integer.parseInt(args[0]) && (one>= Integer.parseInt(args[0])))
{
System.out.println("Both numbers are greater than ceiling. The sum is " + sum);
}
one = one + two;
two = one + two;
}

System.out.println ("The sum is " + sum);
}

**Difference Between Sum of Squares and Square of Sums **

This was my second java script, and it’s again pretty simple. It takes in a user specified number, and calculates the sum of the squares up to this number, as well as the square of the sum of this number, and then returns the difference.

/*
* This script calculates the sum of the squares up to a certain number, n
* and the square of the sums, and finds the difference.
*/

package sumsquarediff;

/**
*
* @author Vanessa
*/
public class Main {

public static void main(String[] args) {
int ceiling = Integer.parseInt(args[0]);
int sum = 0;
int sumsquares = 0;
int i;

//Find the sum of the numbers, and square, and find the squares, and sum

for(i=1;<=ceiling;i++){
sum = sum + i;
sumsquares = sumsquares + (i*i);
}

int squaresum = (sum * sum);
System.out.println("The sum of the numbers squared is: " + squaresum);
System.out.println("The sum of the squares is " + sumsquares);

int difference = squaresum - sumsquares;
System.out.println("The difference is " + difference);
}
}

Largest Palindrome Produced by Multiplying Two Three Digit Numbers

This script finds the largest palindrome produced by multiplying two three digit numbers. Since it was my third java script, I decided to step it up and include a sub-function, called “reverse,” which works recursively to return a reversed string. It works by cycling through potential answers, x and y, starting at 999 for each, and counting down to 1. For each potential pair, it multiples the two numbers, converts this value into a string, reverses it, checks to see if the reversed string is equal to the original number (as a string), and if so, returns the palindrome, which is the answer that we are looking for.

/*
* PaliProd3Nums:  This script finds the largest palindrome
* that is a product of two three digit numbers
*/

package paliprod3nums;

public class Main {

public static void main(String[] args) {
int x;
int y;
int product;
String prodstring;
String reversed;
int largest = 0;
int holder;

for(x=999;x>100;x--)
{
for(y=999;y>100;y--)
{
product = x*y;
prodstring = Integer.toString(product);
reversed = reverse(prodstring);
if (prodstring.equals(reversed))
{
//Tell the user about the possible answer
System.out.println("A possibility is " + x + " times " + y);
System.out.println("which is equal to "; + prodstring);

//Check to see if this answer is larger than our current largest
holder = Integer.parseInt(prodstring);
if (holder > largest)
{
largest = holder;
}
}
}
}
System.out.println("The answer is " + largest);
}
public static String reverse(String input)
{
String holder = null;
String ending = null;

if (input.length() == 1)
{
return input;
}
else
ending = input.substring(input.length()-1,input.length());
String rest = input.substring(0, input.length() -1);
holder = ending + reverse(rest);
return holder;
}
}

Sum Multiples of 3 and 5 Below a Specified Ceiling

This script is pretty dry – it takes in a user input “ceiling,” and then iterates through a loop first checking if the number is a multiple of 3, and then 5. When it finds a number that is either a multiple of 3 or 5, it adds this number to a total sum, and at the end, returns this sum to the user.

function sum_35multiples(ceiling)

%--------------------------------------------------------------------------
% This function takes a number ceiling, finds all the multiples of three
% and five below this ceiling, and then sums them.
%
% Vanessa Sochat     October 23, 2010
%--------------------------------------------------------------------------

% Go through numbers 1 through ceiling, check if multiple of 3 or 5, and if
% it is, add it to our sum.  If it's a multiple of 3, then we skip the
% multiple of 5 check.  If it's not a multiple of 3, we check for 5, so no
% numbers get counted twice.
sum=0;

for i=1:ceiling-1
if mod(i,3)==0
sum = sum + i;
elseif mod(i,5)==0
sum = sum + i;
end
end

fprintf('\n%s%d\n','The sum is ',sum);
end

Format Subject IDs

This script was created to format a list of subject IDs to be used in a python script to batch run an analysis. Currently, the user has to copy paste the IDs from excel into a text document, and then from the text document into nedit, and then spend an inordinate amount of time formatting each ID to be surrounded by parenthesis, and separated by commas. I’ve been the primary runner of batch analysis, which is likely why no one has complained about the process. However, if/when that changes, this will be incredibly useful. With hundreds of subjects, many datasets, and many different types of analysis, doing this manually just… shouldn’t happen! I’m unsure why I didn’t whip this up sooner! /bonk.

So, while this script is incredibly simple, it is incredibly valuable for the lab, so I thought I’d share it.

How Does it Work?

This script takes in a text file, either specified at command line (if you run Format_ID(‘mytextfile.txt’)) or input via a GUI (if you run Format_ID() with no arguments). The text file should be a list of subject IDs, one per line, to be used in the python script. These will be formatted to be surrounded by quotes, and separated by commas. For example, if the text file has the following IDs:

12345_11111

12345_22222

12345_33333

the output will be:

subnums = [“12345_11111″,”12345_22222″,”12345_33333”]

Why a Text File?

I can imagine that some users might prefer a .csv file for input, or perhaps just an excel, however I decided to do this format because I always copy paste my subject IDs from the organizational file into a text file, to be used as reference for various stages of the analysis. So you could say that my reasons were slightly selfish, however I’d be happy to modify the script if anyone has a deep penchant for excel or .csv files.! I, however, will always have a strong preference for good old notepad and wordpad. :O)

You can view the script here

Find all sets in the powerset whose members sum to the largest

The last challenge in the Grelin Code Challenge asks to create a script that finds all subsets of an array where the largest number is the sum of the remaining numbers. For example, if we input (1, 2, 3, 4, 6) we should get the following subsets:

• 1 + 2 = 3
• 1 + 3 = 4
• 2 + 4 = 6
• 1 + 2 + 3 = 6

and the count is 4. The set of numbers that is provided to use is:

• (3,4,9,14,15,19,28,37,47,50,54,56,59,61,70,73,78,81,92,95,97,99)

For this question I couldn’t immediately sit down and start working. I knew that I needed to do some reading up on sets. I first refreshed my memory about powersets – a powerset of a set is all the possible subsets. So, for example, if we have a set (x,y,z), the entire powerset would include: (x),(y),(z),(xy),(xz),(yz),(xyz), plus the empty set (). And very cool – there is an equation to calculate the powerset! You can get the number of subsets in the powerset by using the equation P(s) = 2 to the n, where n is the number of elements in the set. So given that the set provided by the challenge has 22 members, this means that there are 2 to the 22 possible subsets, or 4,194,304 possible answers. Holy cow! I definitely need a script to figure this out!

**So what I decided to do: **

1. I first check for a set that is empty or only has one member. There are no solutions in these cases. If we have more than 1 member, then I cycle through a loop for the number of iterations equal to the number of possible sets in the powerset (2 to the n, where n is the number of possible subsets). We subtract one because the loop starts at 0.

for i = ( 0:(2^numel(theSet))-1 )

3) We then want to find the indicies for each possible set. This means that, for each value of i, we will generate an array of 0’s and 1’s, with a one at each spot where a member is located.

indicies = logical(bitget( i,(1:numel(theSet)) ));

4) We can then use these indicies to get the actual numbers from the set, and put the set into a variable.

set_holder = {theSet(indicies)};

5) We then look at only the subsets that have greater than one member, find the largest value, and compare this value to the sum of the remaining elements. This is where the script would break if the input set wasn’t ordered from least to greatest, because it assumes that the last element in the set is the largest. If the sum matches the largest member, we add it to our candidates variable (p is returned, which feeds into “candidates”), and add one to the count!

% If the set is greater than one element
if length(set_holder{1}) > 1
% Find the index of the largest number, get the number
% (Assuming set organized least --> greatest)
largest_member = find(indicies,1,'last');
largest_member = theSet{largest_member};

total_sum=0;
% Get the sum of the remaining numbers
set_sum = find(indicies(1:length(indicies)),'1');
for z = 1:length(set_sum)-1
total_sum = total_sum + theSet{set_sum(z)};
end

% If the sum of the elements matches the largest number,
% add it to our output set:

if (total_sum == largest_member)
p{candidate_count+1} = {theSet{indicies}};
candidate_count = candidate_count +1;
end
end

6) Lastly, if at the end of this process we haven’t found any answers, we return p as an empty set, so the script notifies the user that there are no sets that fit the criteria.

Things that I would change:

I will post the script that I used to find the answer below, but I also want to talk about some tweaks that I would make to it to make it slightly more efficient. Efficiency and speed isn’t something that I’ve been overly concerned about, but I realize that if I want to get better, I need to start taking it into consideration. So, if I were to add to this script, I would do the following:

• Possibly pre-allocate memory for the variable that holds the answers. I don’t have a great understanding of this stuff, but I think that if a variable changes size on every iteration, it always has to find a new spot in memory. But on the other hand, if I allocate memory based on the largest possible subset, that would be a waste for the smaller subsets. So it’s a question of whether it’s better to take up more memory off the bat, or leave it to allocate on the fly.
• Recursion! There is likely a clever way to do this recursively, and of course that would be even more efficient.
• The challenge just needed the count of the number of subsets that match the criteria, and I decided to save them all into a variable, and then cycle through them to not only return a count, but present the largest to the user. This obviously wasn’t required, and isn’t necessary, and instead I could just present the “candidate_count” variable to the user at the end.
• I think that getting an ordered list was nice, but it is potentially not always going to be the case. So you would want the script to first order the numbers from least to greatest, and then perform the algorithm. That way, we can be assured that the last member of the script is always the greatest. As is, if this script was fed an out of order array, it would still assume the last number to be the greatest, and not function appropriately.

Here is the full script, a la Matlab. How I love Matlab! Keep in mind that this display mucks up greater than and less than symbols.

function powerset(Set)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script takes in a set of numbers, and finds the largest subset
% for which the largest number is the sum of the remaining numbers, and
% returns a count of and list of all the subsets that this criteria applies
% to, as well as the largest subset.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% First we want to cycle through all the subsets of the input set, and save the
% ones that have the largest number be the sum of all other numbers.  Since,
% according to wikipedia, the largest number of subsets for any set is 2 to
% the n, we can go through a loop of that size and create indicies for every
% combination of subsets, and then add subsets to an array that are
% candidates.  We will use the function powerset_review to do this!

candidates = powerset_review(Set);

% Once we have our candidates, we want to find the largest of the set, and
% return this as the answer.

% If we have zero candidates:
if size(candidates) == 0;
fprintf('%s\n','No subsets found for which the largest number is equal to the sum of the remaining numbers')
return
end

% If we have one candidate, return as the answer.
if size(candidates) == 1
fprintf('%s\n','The following subset is the only set in the powerset for which the largest number is equal to the sum of the remaining numbers:')
candidates{1}
return
end

% For all others, when we have more than one candidate, we keep track of the
% largest subset that the pattern applies to, and in the case of finding an
% equivalent length, we keep the first one found:

if length(candidates)>1
fprintf('\n%s\n','All sets for which the largest number is equal to the sum of the remaining numbers include: ')
candidates{1,:}
current_largest = candidates{1};

for k=2:length(candidates)
compare_set=candidates{k};
if (length(compare_set) && length(current_largest))
current_largest = compare_set;
end
end
fprintf('%s%d\n','The number of total candidate subsets is ',length(candidates));
fprintf('%s\n','The following subset is the first largest set found in the powerset for which the largest number is equal to the sum of the remaining numbers:')
current_largest
return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function returns subsets from the powerset that have the largest
% number = the sum of the rest of the numbers!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function p = powerset_review(theSet)
candidate_count=0;

% Deal with empty and one member sets
if isempty(theSet)
fprintf('%s\n','You have provided an empty set.  There is no solution!');
p={};
return
end

if length(theSet)==1
fprintf('%s\n','You have provided a set with only one member.  There is no solution!');
p={};
return
end

% Generate all numbers from 0 to 2^(num elements of the set)-1
for i = ( 0:(2^numel(theSet))-1 )

% Convert i into binary, convert each digit in binary to a boolean
% and store that array of booleans
indicies = logical(bitget( i,(1:numel(theSet)) ));

% Use the array of booleans to extract the members of the original
% set, and check the set to see if the largest number is equal to
% the sum of the additional numbers.  If yes, store the set
% containing these members in the powerset.  This algorithm
% assumes that the set is sorted from least --> greatest.
set_holder = {theSet(indicies)};

% If the set is greater than one element
if length(set_holder{1}) > 1
% Find the index of the largest number, get the number
% (Assuming set organized least --> greatest)
largest_member = find(indicies,1,'last');
largest_member = theSet{largest_member};

total_sum=0;
% Get the sum of the remaining numbers
set_sum = find(indicies(1:length(indicies)),'1');
for z = 1:length(set_sum)-1
total_sum = total_sum + theSet{set_sum(z)};
end

% If the sum of the elements matches the largest number,
% add it to our output set:

if (total_sum == largest_member)
p{candidate_count+1} = {theSet{indicies}};
candidate_count = candidate_count +1;
end
end
end
% If we go through the set and find nothing, be sure to return an
% empty variable
if candidate_count == 0;
p={};
return
end
end
end

and here is a snapshot of finding the answer!

Find all the unique prime factors of a number, and sum them

Here is part two of the coding challenge: to write a script that takes an input number, and finds the sum of all of the unique prime factors of that number. My script is by no means elegant, but it logically goes through a series of steps to get the correct answer, probably in a similar way to how I might solve the problem. It first creates a structural array to hold unique prime factors, and first checks for 0 and 1, and then if the number is divisible by two.

• If we are dealing with 0 or 1 or a negative number. a very prompt is returned!

• If the number is not divisible by 2, we move on to the subfunction “filter_primes” that checks to see if the input number is divisible by the numbers 1 through the (input number / 2 +1). If we find a number that works, we check to see if it’s prime. If it’s prime, then we add it to our array. If not, we move on.

• If the number IS divisible by 2, we just add 2 to our list of prime factors, and then we create a new number, our first number divided by two. We check to see if that number is prime, because if it is, then our answer is simply 2 and that number, and we are done. If the new number isn’t prime, then we go right into “filter_primes.”

At the end, we present the user with the array of unique prime factors, and the sum.

 
function prime_divisors(number) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This function takes in a number, and returns the prime divisors of that % number, and also a sum of those divisors 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% Place the original number in a variable so we can reference it later
original_number = number;
divisor_count = 0;

% First check if we have been given 0 or 1.
if (number == 0 ) || (number ==1 )
    fprintf('%d%s\n',number,' does not have any prime factors.');
    return
end

if (number < 0 )
    fprintf('Please enter a positive number.\n');
    return
end
    
% Check if the number is even (divisible by 2) If it is, divide by
% 2, and save 2 as the first divisor.  If not, use "filter_primes" to 
% find all numbers it is divisible by, and find the list of prime factors.
if mod(number,2)==0
    number = number/2;
    divisors{divisor_count+1}=2;
    divisor_count = divisor_count+1;
    
    % After we divide by 2 to get the new number, we need to check if it's
    % prime.  If it is, then we have finished and have our completed list.
    % We don't want to waste time checking for other potential numbers,
    % because this number = 2 * a prime
    if (isprime(number))
        divisors{divisor_count+1} = number;
        divisor_count = divisor_count+1;
        fprintf('%s%d%s\n','The number ',number,' is prime, and we have finished searching for primes.  Adding to list!')
        fprintf('%s\n','The list of primes is:')
        
        % Print out the list of prime divisors, as well as a sum.
        report_answer(divisors);
    else
        filter_primes(number);
    end
else
    filter_primes(number);
end
    
    
%--------------------------------------------------------------------------
% Function isprime
%
% This function checks if an input number is prime by dividing by all
% possible factors through the input number / 2.  In the case that there 
% is always a remainder for all these divisions, we know that the number
% has to be prime, meaning that it is only divisible by 1 and itself.
%--------------------------------------------------------------------------
function f = isprime(number_two)

    % if our number is even, we don't continue
    if mod(number_two,2)==0
        fprintf('%s%d%s%d\n','The number ',number_two,' is even so it cannot be prime');
        f = 0;
        return;
    else
        for i = 3:(number_two/2)
            if mod(number_two,i) == 0
                fprintf('%d%s%d%s\n',number_two,' is divisible by ',i,' and therefore cannot be prime!');
                f = 0;
                return;
            end
        end
        f = 1;
        fprintf('%d%s\n',number_two,' is not divisible by anything and is therefore prime!');
        return
    end
end
%--------------------------------------------------------------------------
% Function filter primes
%
% This function cycles through a number and finds all the numbers it is
% divisible by.  When it finds a number that it is divisible by, if the
% number is prime, it adds it to the list we are saving.
%--------------------------------------------------------------------------
function h = filter_primes(number_three)

    fprintf('%s%d%s\n','The number ',number_three,' is not even, checking divisibility')
    for j = 3:((number_three/2)+1)
        if mod(number_three,j) == 0
            fprintf('%d%s%d%s%d%s\n',number_three,' is divisible by ',j, ' checking if ',j,' is prime...');
            if (isprime(j)==1);
                divisors{divisor_count+1}=j;
                divisor_count=divisor_count+1;
            end
        end
    end
     report_answer(divisors);   
end
%--------------------------------------------------------------------------
% Function report_answer
%--------------------------------------------------------------------------
function r = report_answer(divs)

fprintf('\n%s%d%s\n','The list of prime divisors for ',original_number,' is:')
divs

sum=0;
for z= 1:divisor_count
    sum = sum + divs{z};
end
    
fprintf('%s%d\n','The sum of this list is: ',sum)
return
end
end

Here is a visual. I could remove a lot of the print lines (minus the answer) – these were just for debugging purposes!

Find the smallest prime greater than an input number

that is also in the Fibonacci Sequence!

For this one I left in all my print lines that were used for debugging / testing purposes. This script takes in a number, and finds the smallest prime that is greater than that number that is also part of the Fibonacci Sequence. It starts with the numbers 1 and 1, of course (the first two in the Fibonacci Sequence) and moves upward in the sequence until it hits a member of the sequence that is above the input number. It stops at the two members of the sequence before this number, and then adds these two numbers, tests if the outcome is prime, and if so, returns the value as an answer. If not, it calculates the next number in the sequence, and continues testing.

• The main function is called prime_fibonacci() with the number going in as input. There are two subfunctions one_up, which moves us up one order in the sequence to the next two numbers (so from 1, 1 to 1 2, for example) and

• is_prime: which does exactly what you would think – it first tests if the number is divisible by two, and if so, returns false (0). If not, it takes in an input number and tests for primeness by dividing by all possible factors – 1 through the (number + 1) /2. The modulus function is used for this test. In the case that there is no remainder, this means that the number is not prime, and we return false. In the case that there is always a remainder, then hooray! The number is prime, and we return it as the answer.

function prime_fibonacci(number)

% Find 2 fibonacci numbers below the number given by the user.  Start with
% n_one as 1, and n_two as 1, and slowly work upwards until we are just
% below the user specified number.
n_one = 1;
n_two = 1;

% n_one is always larger, we check it first.  At the end of this loop we
% have the two numbers in the Fibonacci sequence directly before our input
% number
fprintf('%s%d%s%d\n','Number one is ',n_one,' and Number two is ',n_two);

while n_one < number && n_two < number
if n_one + n_two <= number
[n_one,n_two] = up_one(n_one,n_two);

fprintf('%s\n','The summed numbers are still less than the input number, getting next set.');
fprintf('%s%d%s%d\n','Number one is now ',n_one,' and Number two is ',n_two);
else
fprintf('%s\n','The summed numbers would be equal or greater than the input number, stopping');
break;
end
end

% We add the two numbers to get the next number in the Fibonacci sequence,
% which we know is larger than our number after the tests above.  We now
% need to check if this number is a prime.  If yes, then we've found our
% answer.  If not, then we need to look at the next number.

decision = 0;
while decision == 0;
decision = isprime(n_one + n_two);
if decision == 0
fprintf('Creating a new set of numbers...')
fprintf('%s%d%s%d\n','New Numbers: Number one (',n_one, ') and Number two: (',n_two, ')');
[n_one,n_two] = up_one(n_one,n_two);
else
answer = n_two + n_one;
fprintf('%s%d%s%d%s%d\n','The sum of Number one (',n_one, ') and Number two (',n_two, ') is a prime number! ',answer);
fprintf('%s%d\n','The answer is ',answer);
end
end

%--------------------------------------------------------------------------
% Function isprime
%
% This function checks if an input number is prime by dividing by all
% possible factors through the input number / 2.  In the case that there
% is always a remainder for all these divisions, we know that the number
% has to be prime, meaning that it is only divisible by 1 and itself.
%--------------------------------------------------------------------------
function f = isprime(number_two)

% if our number is even, we don't continue
if mod(number_two,2)==0
fprintf('%s%d%s%d\n','The number ',number_two,' is even so it cannot be prime');
f = 0;
return;
else
fprintf('%s%d%s\n','The number ',number_two,' is not even, checking divisibility')
for i = 3:(number_two/2)
fprintf('%s%d\n','Checking divisibility for ',i)
if mod(number_two,i) == 0
fprintf('%d%s\n',number_two,' is divisible by ',i, ' and therefore cannot be prime!');
f = 0;
return;
end
end
f = 1;
fprintf('%d%s\n',number_two,' is not divisible by anything and is therefore prime!');
return

end
end

%--------------------------------------------------------------------------
% Function upone
%
% This function creates the next two fibonacci numbers.  n_one is always
% treated as the larger of the two.
%--------------------------------------------------------------------------
function [numone,numtwo] = up_one(one, two)

numone = one + two;
numtwo = one;
return

end

end

Find the Longest Reversed Substring in a String

This function takes in a string, “stringy,” and finds the longest reversed substring (palindrome) in the larger string. So, racecar would be an example because it is the same thing backwards and forwards!

I basically did it by identifying all the possible centers in the string (indicated by the pattern xyx, and then I went through those, and looked at the character directly to the left, and directly to the right, and checked for a match. If it matched, I then went one letter farther out, until they no longer matched, and the longest string found for each is recorded in a structural array with two fields: the answer, and the length. When I was finished looking through all of my possible centers, I print out the one with the longest length, and that is the answer!

The string given for the challenge was:

'FourscoreandsevenyearsagoourfaathersbroughtforthonthiscontainentanewnationconceivedinzLibertyandded
icatedtothepropositionthatallmenarecreatedequalNowweareengagedinagreahtcivilwartestingwhetherthatnapti
onoranynartionsoconceivedandsodedicatedcanlongendureWeareqmetonagreatbattlefiemldoftzhatwarWehave
cometodedicpateaportionofthatfieldasafinalrestingplaceforthosewhoheregavetheirlivesthatthatnationmightliv
eItisaltogetherfangandproperthatweshoulddothisButinalargersensewecannotdedicatewecannotconsecratew
ecannothallowthisgroundThebravelmenlivinganddeadwhostruggledherehaveconsecrateditfaraboveourpoorpon
wertoaddordetractTgheworldadswfilllittlenotlenorlongrememberwhatwesayherebutitcanneverforgetwhatthey
didhereItisforusthelivingrathertobededicatedheretotheulnfinishedworkwhichtheywhofoughtherehavethusfars
onoblyadvancedItisratherforustobeherededicatedtothegreattdafskremainingbeforeusthatfromthesehonored
deadwetakeincreaseddevotiontothatcauseforwhichtheygavethelastpfullmeasureofdevotionthatweherehighly
resolvethatthesedeadshallnothavediedinvainthatthisnationunsderGodshallhaveanewbirthoffreedomandthatg
overnmentofthepeoplebythepeopleforthepeopleshallnotperishfromtheearth'

and my scripty-doo got it right on the first try! :O)

function reverse_substring(stringy)

string_length = length(stringy);

count = 1;

for i = 1:string_length-2
if strcmp(stringy(i),stringy(i+2))
center{count}= i+1;
count = count+1;
end
end

for i = 1:size(center,2)
if string_length/2 >= center{i}

extent = center{i};
else extent = string_length - center{i};
end

j=1;
while j < extent-1
if strcmp(stringy(center{i}+j),stringy(center{i}-j))
answer = stringy(center{i}-j:center{i}+j);
j = j+1;
else break
end
end
CONTENDERS(i) = struct('Answer',answer,'Length',(2*j)+1);
end

longest = CONTENDERS(1).Length;
for i = 1:size(CONTENDERS,2)
if CONTENDERS(i).Length > longest

longest = CONTENDERS(i).Answer;
end
end

longest

end

Automated Reminder Emails and Text Messages with Matlab and Batch

I’ve worked about 35 hours in the past four days putting together an automated reminder system for our lab. I basically have a batch script working with two matlab scripts, GMail calendars, and various data files to send automatic email and SMS reminders to participants in our study. The system isn’t yet implemented, but I’ve tested it and have everything in place to launch.

Where is it run from?

All scripts are located in a secure script directory in an Appointments subdirectory. The main batch script that coordinates the entire process is executed via a scheduled task, which would be set to run nightly in the wee hours of the morning. The requirements of the computer running the script are 1) access to the Appointment folder and 2) an installation of matlab.**

How does it work?

This Appointments Folder contains scripts and files necessary for sending out automatic, nightly reminder emails, text reminders, and email confirmations for imaging and computer battery data. This documentation will explain the contents of the folder, as well as the workings of the various scripts.

Overview of Files

rem_email.bat This batch script is responsible for downloading necessary appointment data from the private server, and running the matlab scripts to send out reminder and confirmation notifications. The script would be run nightly to ensure that new files have been downloaded and we have the most up to date information. It basically maps a drive to download new data, runs the matlab scripts to send out reminders, deletes temporary files produced by the scripts, and then closes the mapped drive. The script is as follows (with all paths and such replaced with text):

rem_email.m

schedule.csv: Is scheduling information saved on a private server which we connect to and download for the scripts.

Takes in the following arguments: rem_email(days,cal,sendtxt,sendemail,contactemail)

• days: This is the number of days in advance that the reminder emails will be sent. (1 = 24 hrs, 2 = 48 hrs, etc)

• calendars: This is the calendar name that the reminders will be sent to, for each of the times indicated. This can currently be “Computer” or “Imaging.”

• sendtxt: ‘yes’ indicates we want to send a text message, no means no!

• sendemail: ‘yes’ indicates we want to send an email, no means no!

• contactemail: (optional) An email address that you want to be notified in the case of error

This script depends on the imaging and computer calendars to be formatted with the participant first and last name in the Title (Summary) field, and the appointment information copied from the clinical calendar into the description field the following format:

Name: Jane Doe
Phone: (999) 999-9999
E-mail: jane.doe@duke.edu

The script works as follows:

1. Downloads the latest computer or imaging calendar based on a private ical address, and saves the data as either ‘imaging.txt’ or ‘computer.txt.’
2. Reads schedule.csv (just downloaded) into matlab. In the case that this file doesn’t exist (due to a download error, etc) it records this into the error log, and contacts the ‘contactemail’ about the error.While reading in information from the calendar file, the script places the participant name, and appointment date and time into a structural array ONLY for appointments in the future.
3. All of this data is read in from particular locations in the calendar file, and formatted between numbers and strings so we can both do calculations (numbers) and print the variables (strings) into text in reminder emails, etc.
4. Important Notes about the formatting of the calendar: In the case that the script cannot find the person’s name as the first field under “Description” – it extracts it from the “SUMMARY” field – which is like the Title of the event on the calendar. If both these locations are missing the participant name, the script will exit with error
5. The script next uses the fields from the schedule.csv file as a lookup table, and finds participants by first and last name to match with individuals with future appointments only. Currently, we have no third check, so if two people have the exact same name, the script will find the one that appears earliest in the file. This hasn’t been an issue thus far, but if it becomes an issue I will troubleshoot a solution.
6. As the script matches participants in the structure (read in from the calendar file) with additional information read from the schedule.csv, it updates the structure, so at the end of this process we have a complete structural array with information for all subjects with future imaging or computer appointments. Since these are scheduled always within a week of the clinical appointment and never greater than a week apart, this structure should never get greater than perhaps 30. Even if it did, it wouldn’t be an issue.
7. Once we have this complete structure, we then need to figure out who to contact.
8. The script looks at the “days” variable, specified by the user, which is the number of days in advance to send the reminder for. It converts the date of the present date into a number, adds the “days” variable to that, and creates a “compare_date” variable.
9. We then cycle through the structural array of participants, and convert each person’s appointment date to a number, and then compare this number to the comparedate. In the case that we have a match (meaning that the participants appointment is indeed X (days) away, then we send the email.
10. If the person has indicated that it’s OK to text message them, and the script is set to sent text reminders (determined by the sendtxt variable) then we also use the script send_text_messages to send the text reminders.
11. Depending on whether we are doing the imaging or computer calendar, the email message varies in telling the participant the appointment type (“Imaging” or “Computer”) and the location to meet the experimenter.
12. We then record all participants contacted to an output log, and exit matlab to prepare for the next script.
13. The rem_email.bat file can send out different combinations of reminder emails and texts simply by running this script multiple times with various configurations, as is done if you look at the batch code above.

confirmation_email.m

Takes in the following arguments: confirmation_email(contactemail)

This script works in the same basic manner as rem_email.m, except that it is hard coded to download the imaging and computer calendar, and then send out only reminder emails to subjects with appointments that were just created in the past 24 hours (new appointments that should be confirmed). It figures out this detail by reading in the CREATED field from the calendar text files, which contains the date and time of when the event was created on the calendar. Output goes to a second output log, and in case of error, the contactemail is notified.

imaging.txt: Is a temporary file created in an ical format (saved as a text for easy readability) that contains information from our imaging calendar on Gmail. This file is created by the script (if it doesn’t exist) and deleted at the end of the batch job.**
**

**computer.txt: ** Is the equivalent temporary file, but for the computer battery calendar.
various LOG.txt files: Are output logs and error logs of reminder emails and texts sent. To be updated upon each script run.

These are new scripts, of course, and should be checked regularly for successful runs, and in the case of error, troubleshooted! In the long run I am hoping these will provide an easy and reliable toolset for sending reminder emails and text messages to participants, and doing away with missed appointments!

Convert Military Time to User Friendly String

I want to share this, because it could be useful in many contexts. I am working on a script that downloads information from a Google Calendar, and automatically sends reminder emails based on a user specified time. Given that the times are military, you don’t exactly want the email to say “You have an imaging appointment at 1730,” so I wrote this quick little portion of my script to, after pulling the time from the text file, convert it to a user friendly string, like “5:30 pm.”

This script starts when it knows that it’s on the line that contains the appointment time (the line is stored as the “currentline” variable, and for this line in the file exported from the gmail calendar, the actual time always starts and ends at the same character location.

% read in military time
time = currentline(18:21);
% Convert military time to standard time
if (str2double(time)>1259)
time = str2double(time)-1200;
period = 'pm';
else
period = 'am';
end
% Format time by separating the last two characters
% from the first, and sticking them together with the :
time_end = time(length(time)-1:length(time));
time_beg = regexprep(time, time_end, '','once');
% Put it all together into a user friendly format for printing
time = [ time_beg ':' time_end ' ' period ];

And now I can save it into a structural array of subjects, and use it when I create the text for my email. Hooray! Also note that the “>” should be a greater than symbol, I’m not sure why it’s stubbornly coming out like that!

Get rid of carriage returns using tr

Kristin recently had a script that, for some mysterious reason, was in Windows format! (the error output obviously indicated the presence of carriage returns!). In my troubleshooting, this was the command that fixed it!

• tr -d ‘\15\32’ fixed_script.sh

Plink_View (Plink_View.m) Alpha Release!

When doing any sort of genetic analysis, we use a program called Plink to read in subject data from a .ped/.map file, and help us pull the SNPs that we are interested in. Other variables in Plink include individual and family IDs and the allele combo that each participant has for each SNP. It was a pretty cumbersome process to take the output produced by this program (a strangely formatted .list file) and wrestle with it to get it into excel columns, and then come up with some custom (likely error prone) process to match the many subject family and individual identifiers with whatever identifier is used in the study. This script makes that process, which might have taken someone an hour, complete in under a minute!

This script reads in a .list file created with the software program “Plink” and allows the user to select SNPs that he/she is interested in, as well as a subject ID lookup table, and a gene reference table. The script matches the study IDs with the correct participants, as well as appends the gene identifier to the name of each SNP. The output is a .csv file for easy import into Excel, SPSS, etc. The name of the output file is the name of the list file appended with the date.

Running the Script

1. Save the script to your directory of choice, and prepare a .list file, and optionally, a subject lookup table and a genetic reference file.
2. Either add the location of the script as a path (addpath(genpath(‘mypath’)), or navigate to where the script is to run it.
3. To run, simply type “Plink_View” on the MATLAB command line. You can either run with no input arguments, or specify your listfile and lookup table as arguments 1 and 2, respectively. Here are your choices:

• Plink_View – will be prompted for listfile, lookup, and ref

• Plink_View(listfile) – will not have gene names and other_IDs

• Plink_View(listfile,lookup) – will not have gene names

• Plink_View(listfile,lookup,ref) – will have gene names and other_IDs

Variables

• ** listfile**: Must be a .list file generated by plink, with each row listing the gene identifier, the SNP name, the allele combo, and then groups of family IDs followed by individual IDs. The script expects various allele combos for each SNP (11, 12, 22, 00) and family IDs to be followed by the individual ID.

• **lookup: ** Must be an excel file (.xls/.xlsx) with two columns – it is OPTIONAL. The individual identifier should be in the first column, and the “other_ID” in the second column. The user is allowed to select if he/she would like to lookup by individual ID, family ID, or Family_Individual ID (the concatenation of the two).

• ref: is an rs_gene_reference excel (.xls/.xlsx) document that is OPTIONAL. It should be formatted so that the first column has the SNP name, and the second column has the gene symbol. If the user specifies a gene reference file, then the gene symbol is appended to the SNP name in the output table. If not, the output table will only have the SNP name.

How does it work?

**1. The script starts by evaluating how much information the user has provided. In the case of zero arguments, it prompts the user for all of the input variables (giving the choice to add the lookup table and gene reference table, since those two are optional!) In the case that input arguments are provided on the command line, it goes right into step 2.

1. The script reads the second column of the list file, which is always the long list of SNPs. Regardless of the number of repeats of each SNP (based on the number of allele combos) the script presents the user with a popup GUI with a complete list of the SNPs. The user is allowed to select the SNPs from the list that he/she wants included in his/her output file.

1. The script then loads the genetic reference table, and matches each chosen SNP with its gene. These names are saved into another variable for use in the creation of the output file. They will be the headers for each SNP. For example, if a SNP is called rs16 and it is from the gene “THSD7A”, the column header that includes each subject’s allele will have the title “rs16_THSD7A.” In the case that the user does not use a gene reference table, the name will simply be “rs5522.

2. The script next delves into the task of reading the data. This is reliant on Plink producing a consistently formatted output .list file. Here is a sample output with two SNPs (shortened of course!):

7 rs16 11 7612 1 6931 1 3658 1 4346 1 9744
7 rs16 12 941 1 8955 1 9705 1 7220 1 9232 1 5650 1 698 1
7 rs16 22 5966 1 9412 1 9519 1 5922 1 2772 1 7129 1 5537 1
7 rs16 00
7 rs8 11 7680 1 7468 1 3862 1 4992 1
7 rs8 12 8955 1 9519 1 5922 1 4346 1 8796 1 8905 1 8702 1
7 rs8 22 941 1 7612 1 5966 1 6931 1 9412 1 9705 1 3658 1 7220 1 \ 7 rs8 00 6904 1 7312 1 7393 1

Each line follows the following format:

(gene) (SNP) (allele) (family ID 1) (individual ID 1) (family ID 2) (Individual ID 2)… etc

1. The script reads in one line, one element at a time, skipping over the (gene). This is done by locating the next empty space, and then referencing the line from the start to that space -1. When an element has been read and saved to a temporary variable, it is cut from the beginning of the line.

2. The SNP name is added only in the case that the user indicated wanting it. A DATA structure is created that holds the various fields of each person. In the case that we have a .list file with multiple SNPs, it is very likely that an individual will appear twice. So each person is read in with their family and individual ID, and then searched for in the structure. In the case that they already exist, the new SNP information is added. In the case that they don’t exist, a new subject is added with the new family and individual ID.

3. Once all subjects and their SNPs have been added, the script loads the study ID lookup table, and goes through each subject saved in the DATA structure, searching for a matching ID in the lookup table. In the case that the ID is found, it is placed into the “Study_ID” field for the subject. In the case that it is not found, this field will remain blank.

4. The last step is printing all of the data to a .csv file with the proper headers and comma placement. The script gives the .csv file the same name as the list file provided by the user appended with the date.

View the Script

"Send Brains" (send_brains.m) Alpha Release!

It’s a weekly duty in our lab to preprocess data and send out a brain image (with a lab logo) for each subject that has participated in our study. It used to be one of my weekly duties, and I handed my scripts over to an undergrad to take charge of. However, it stayed with me how tedious the entire process was, and I added it to my list to create a better way, when I found the time!

I’m happy to announce the release of send_brains.m – a matlab script that I wrote to automatically create AND send brain images. While this isn’t like the coverage checking script in that it would be easy to run “as is” for another lab, with some modifications it could easily be used. It has a couple really cool features, including using an SPM function called “slover” to prepare axial and sagittal views of a high resolution anatomical, saving .jpg and .zip files, and sending an email directly from MATLAB.

What Does it Do?

The script looks at subjects that have been Processed (through the BOLD processing pipeline that I created for our lab), and compares this list with completed brain images in a folder where all brain image output is stored. Subjects that are in the Processed folder but not in the Brain Images folder are added to a list of potential sendees.

How Do I Run It?

• To run the script, make sure that you have its path added. We have a Script folder with many subfolders, so if you have something similar, it’s easiest to add the path with subfolders : addpath(genpath(‘Your/path/here/’));
• To run, simply type brain_images on the MATLAB command line.

How does it work?

1. The script first creates the list of potential sendees – whoever has been newly processed that week and does not have a brain image created yet.
2. The script then presents this list to the user in a graphical window, and the user can select subjects to send brain images to by highlighting the user IDs. (updated 10/12/2010) ( The old method asked the user if he/she would like to delete any participants, and the user had to enters the ID into the GUI, and then the script would present the user with the new list. I wanted to make this easier so I changed it.)
3. Next, the script looks for the processed anatomical in the subject’s anat directory under Processed. This image will normally always have the same name – however in rare cases when the name is different, for whatever reason, the script will not find the highres, and will present the user with a GUI to select the image that he/she would like to use for that participant. It then prepares two slice views of the highres, crops them, puts them together as a zip, moves the zip to the Graphics/Brain Images folder, and cleans up the old files.
4. Once all brain image zips have been created, the script prompts the user if he/she wants to send a brain image for each subject, and asks the user to enter the email. This process is fairly rapid and easy. In the case that the user mistypes an email, the easiest thing to do is delete the output image, and quickly run the script again.
5. As each address in entered, the script sends the email directly from MATLAB.

Instructions

**1. **Run by typing “send_brains” into MATLAB. Select your experiment top directory.

2. Select your brain images output folder:

3. **Select the subjects that you want to send images to from the GUI. In the case that the standard highres is not present for a subject, the script will prompt you to choose another file, or select “cancel” to skip the subject.
**4. **Go have a sandwich. The script will now first create a slices image with a logo overlay for each participant…

**5. **…and then a sagittal slices view.

**6. When everyone’s images have been created and intermediate files cleaned up, then you will be prompted for each subject if you want to send a brain image or not:

**7. **and then lastly you enter the participant’s email. That’s it!

The one con (between this method and the old manual one) is that we are moving from processing on the cluster (parallel) to a local machine (one at a time!). So instead of everyone being processed at once, and the images taking 30 seconds for everyone, the time to send the images is n X 30 seconds. However, the key point is that you don’t have to do anything! With the old method, you would have to log into the cluster, navigate to the script, manually figure out the new IDs to run and input them into the script, chmod u+x and run it, then manually open up each image in FSLview, take a screen shot, copy paste into paint, copy the logo from a URL, paste and resize it, and then save the entire thing as (Subject_ID).jpg. Yep, I’d say that the script is a heck of a lot easier!

The scripts that are needed include
send_brains.m

crop.m (I did not write this script, but it is necessary to crop the images)

If you would like to see sample output, look at:
Zip Send to Participant

Once and Forever

Vibrance and rapture, epiphany and sound
Gummy and splendor, elegance and bound
Constant in motion, mind in a whirl
Responsibility of a woman, heart of a girl

Glowing the night lamp, hard fall the eyes
Gentacular hum of the heater subside
Blanket the darkness, creep out the mind
Uncover that which to day is blind

Enveloped by passion, research abound
No chance to pause, no heart to be found

But near draw the night, a moment of clear
Escape from the heart a memory so dear
Four years past, sunshine and snow
Imprinted in memory, a face does not go

Quiet in weakness, imprinted by heart
Destined to find, and be jostled apart
A chance of a voice, a hint of a scent
Painful to suppress that which is meant

Never again, vows before dream
Comfort in one, spot, bed, and seam
Anew springs the day, no time to feel
Covered up gently, this Achilles heel

Stronger to face, resilient to go
Naive to unwilling, I never shall know
Sunsets in play, day night and sky
That which is strong wrenching to cry

Burgeon my senses, mind, soul and heart
loving you forever, though forever apart

PickAtlas Version 3.0 Released for SPM

PickAtlas, an SPM extension created by the folk over at Wake Forest, has just had the release of Version 3.0. I’ve installed the new version in SPM8, and wanted to give everyone some feedback. The biggest change is with the addition of masks – there are now rodent atlases, monkey atlases, and other small critters that my lab is unlikely to use. However, it’s very cool that we have them, if we ever decided to look at critter brains!

An important change to note is the use if nifti images from analyze. All of the atlas image files are nifti, as well as any mask that you save from the PickAtlas GUI. When I first tried to add the JHU White Matter Atlas (that I had added to my older version) it spat out a big ugly error – which was fixed when I converted the .img/.hdr file into a .nii. So if you want to add any atlases, the procedure is exactly the same, but you must use nifti!

When using it for “Results” – the only slight tweak is user experience:

The only change that will be noticeable to the user from this point on **(after selecting PickAtlas GUI or a mask from file) **is that the mask will be resliced and a completion meter will be displayed during this process. Because mask area is generally smaller than the whole brain, the number of multiple comparisons will be reduced. Thus, results viewed using the PickAtlas Tool will include a small volume correction that will be reflected in the p-values.

I’m not sure if this was there before – but there is a “Generate Table” button that allows you to select one or more analyze images with an ROI, and it spits out a table with the following statistics:

• Size
• Average
• Std.Dev.
• T
• Region
• ROI name
• L/R/B
• Study
• Image
• Max Value
• Max Loc
• Min Value
• Min Loc

These definitely might be useful!

Everything else is the same – it still does talarach/MNI coordinates, allows for dilation and shapes, right vs left selection, all that good stuff!

Download PickAtlas 3.0 Here

Happy PickAtlasing everyone!