text:
# Implementing k-nearest neighbors in python

Last time, we looked into one of the simplest classification algorithms in machine learning called binomial logistic regression. In this post, I am going to cover another common classification algorithm called K Nearest Neighbors, otherwise known as KNN.

To recap, we have mostly discussed regression models such as simple and multivariate linear regression and polynomial regression which are used for predicting a quantity. On the other hand, classification models are used for predicting a category such as yes/no, will buy car/scooter/truck, will turn pink/green/red etc.

KNN is another such classification model and once I am done explaining how it works, you will understand from where it gets its name. KNN is a supervised learning algorithm that follows 5 simple steps to predict a class:

1. Pick a value for k – this is the number of neighbors you will consider for each new datapoint you want to predict a class for.
2. Pick the neighbors – calculate distance from new datapoint to rest of the data points. The distance can be euclidean distance, manhattan distance or something else.
3. Identify categories – Pick the k data points that are closest to new datapoint and then identify which category they belong to.
4. Label – Label the new datapoint category which most k nearest data points belong to.

Key points to remember about KNN algorithm are that:

• KNN doesn’t make any assumptions about initial dataset unlike some other popular machine learning algorithms. This means that KNN can work with more diverse datasets and hence, can be more accurate.
• KNN computes the results on the go. It memorizes the training dataset and then uses that to label the new dataset when it is run.
• KNN doesn’t scale well. It can be very slow for large datasets.

## Exploring the dataset¶

For this example, we are going one of the most popular datasets in machine learning called the Iris dataset. Some of you might have heard of this dataset already. This data was collected by Fisher in 1933. It contains 4 features and a label. The features are:

• sepal length
• sepal width
• petal length
• petal width

and the label is the type of iris plant.

Our task is to use these four features to correctly predict the type of iris plant.

In [1]:
# Let's load the data into python and take a look at it
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

Out[1]:
sepal_length sepal_width petal_length petal_width class
0 5.1 3.5 1.4 0.2 Iris-setosa
1 4.9 3.0 1.4 0.2 Iris-setosa
2 4.7 3.2 1.3 0.2 Iris-setosa
3 4.6 3.1 1.5 0.2 Iris-setosa
4 5.0 3.6 1.4 0.2 Iris-setosa
In [2]:
# Let's modify the class values so that they don't have "Iris-" prefix
dataset['class'] = [x.split('-')[-1] for x in dataset['class']]

Out[2]:
sepal_length sepal_width petal_length petal_width class
0 5.1 3.5 1.4 0.2 setosa
1 4.9 3.0 1.4 0.2 setosa
2 4.7 3.2 1.3 0.2 setosa
3 4.6 3.1 1.5 0.2 setosa
4 5.0 3.6 1.4 0.2 setosa
In [3]:
# Let's get some more information about our dataset
dataset.describe()

Out[3]:
sepal_length sepal_width petal_length petal_width
count 150.000000 150.000000 150.000000 150.000000
mean 5.843333 3.054000 3.758667 1.198667
std 0.828066 0.433594 1.764420 0.763161
min 4.300000 2.000000 1.000000 0.100000
25% 5.100000 2.800000 1.600000 0.300000
50% 5.800000 3.000000 4.350000 1.300000
75% 6.400000 3.300000 5.100000 1.800000
max 7.900000 4.400000 6.900000 2.500000

We can see there are 150 rows. The mean value of sepal_length, sepal_width, petal_length and petal_width is 5.84 cm, 3.05 cm, 3.76 cm and 1.20, respectively.

We don't see information about our label because it's values are not quantitative.

## Preprocessing the dataset¶

Iris is a very standard dataset and doesn't have any missing or null values so we don't need to do any preprocessing.

Only thing we need to do is split the dataset into dependent variable and independent variable so that we can feed it to our machine learning model.

Keep in mind, that life is never this easy and data preprocessing is a very crucial step in the machine learning process and often the most time consuming.

In [4]:
X = dataset.iloc[:, :-1].values
y = dataset.iloc[:, -1].values

## Splitting the dataset into training and testing set¶

KNN is a supervised learning model (see my earlier post on supervised and unsupervised models) which means we need to feed it some data first to train it. And, once it is trained, we need to test it on a different set of data to evaluate it. To do this, we need to split our dataset into training and testing set. Rule of thumb is to assign 80% of the dataset to training and 20% to testing.

In [5]:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2)

At this point, we have four different arrays:

• X_train - independent variable (training set)
• X_test - independent variable (testing set)
• y_train - dependent variable (training set)
• y_test - dependent variable (testing set)

## Building the model¶

This is the step where we actually build our model.

In [6]:
from sklearn.neighbors import KNeighborsClassifier

# Create classifier using K Neighbors Classifier with k=5
classifier = KNeighborsClassifier(n_neighbors=5)

# Training our model
classifier.fit(X_train, y_train)

# Predicting values using our trained model
y_pred = classifier.predict(X_test)

## Evaluating the model¶

We now have our predicted values, y_pred, and our actual values, y_test. Let's look at how well we predicted.

In [7]:
# Generating a confusion matrix
from sklearn.metrics import confusion_matrix
conf_matrix = confusion_matrix(y_test, y_pred)

print(conf_matrix)

[[ 8  0  0]
[ 0 12  2]
[ 0  0  8]]

Let's calculate accuracy score which is the most popular metric for classification models:

In [8]:
from sklearn.metrics import accuracy_score
print(accuracy_score(y_test, y_pred))

0.933333333333

With a 93.3% accuracy score, I would say our model is quite accurate. We can improve the performance more by fine tuning the hyperparameters such as the initial value of k we selected earlier on to build our model. A great way to decide which values to pick for our hyperparameters is known as cross validation which we will cover in a different post later.

Hope you found this post useful. You can download this code from my github.


text:
# Standalone Win32 for v1.84.0 doesn't include pygame

My experiment that used pygame* for presenting sound runs fine in v1.83.04 but fails in v1.84.0 telling me to install an audio library. In case anybody else runs into this problem, I resolved it by copying over the pygame directory from v1.83.04’s C:\Program Files (x86)\PsychoPy2\Lib\site-packages directory into the same directory for v1.84.0.

Is this a packaging error or has pygame support been deprecated?

• pyo isn’t usable for me because I need it for voicekey functionality, in which case PsychoPy will crash with an obscure error if the audio library preference isn’t set to just [‘pygame’]

Ah, sorry about that. I had to recreate the build package for win32 and obviously forgot (I haven’t used pygame in a logn time). I’ll add it next time

1 Like

ok, thanks! Good to know I didn’t simply screw up something in my experiment code

Dear @jon,
It seems Standalone PsychoPy 1.84.2 win32 also doesn’t include pygame!
Could you explain how I can install pygame? I’m quite a newbie with all of this… I use the builder view and am inexperienced in programming…

Thank you!!
Raina

P.S.: The problem I have is a windows error that pythonw.exe has stopped working, see this thread: http://discourse.psychopy.org/t/windows-error-pythonw-exe-has-stopped-working-please-help/1388 If you can also contribute to this original issue, that would also be much appreciated!

And this is where I got the idea that changing the order in the Psychopy preferences to chose pygame over pyo could maybe help…: http://discourse.psychopy.org/t/randomly-python-exe-has-stopped-working-after-windows-update/1232/10 Do you think it might work in this situation?

format: plain


text:
# Practical problem solving¶

So far we have been learning about computing and coding.

The purpose of this chapter is to introduce some practical techniques for tackling data analysis problems.

At this point we go back to looking at the Swiss-Prot FASTA file used in First steps towards automation. You can download this file using the curl command below.

$curl --location --output uniprot_sprot.fasta.gz http://bit.ly/1l6SAKb  The downloaded uniprot_sprot.fasta.gz file is gzipped. You can view the content of this file without permanently unzipping it using the command below. $ gunzip -c uniprot_sprot.fasta.gz

If the curl and gunzip commands above look unfamiliar you may want to review the First steps towards automation chapter.

By briefly inspecting some random FASTA description lines it becomes clear that some species include variants. Below are some sample E. coli extracts from the FASTA file.

...IraP OS=Escherichia coli (strain SE11) GN=iraP PE=3 SV=1
...IraP OS=Escherichia coli (strain SMS-3-5 / SECEC) GN=iraP PE=3 SV=1
...IraP OS=Escherichia coli (strain UTI89 / UPEC) GN=iraP PE=3 SV=1

How many variants are there per species? Do all variants contain the same number of proteins?

To answer these questions we will make use of the Python scripting language and some practical problem solving techniques.

## Stop and think¶

One of the most common mistakes when programming is to start coding before the problem is clearly defined.

What do we want to achieve? Pause for a second and try to answer this question to yourself. It may help trying to explain it out loud.

Now describe all the steps required to achieve this goal. Again it may help speaking out loud to yourself or imagining that you are describing the required steps to a colleague.

Did you find that easy or difficult? Thinking about what one wants to achieve and all the steps required to get there can often be quite overwhelming. When this happens it is easy to think to oneself that the detail will become clear later on and that one should stop wasting time and get on with the coding. Try to avoid doing this. Rather, try to think of a simpler goal.

For example consider the goal of working out how many different species there are in the FASTA file. This goal would require us to:

1. Identify all the description lines in the FASTA file
2. Extract the organism name from the line
3. Ensure that there are no duplicate entries
4. Count the number of entries

If you know the tools and techniques required to do all of the steps above this may not seem overwhelming. However, if you don’t it will. In either case it is worthwhile to treat each step as a goal in itself and work out the sub-steps required to achieve these. By iterating over this process you will either get down to steps sizes that do not feel overwhelming or you will identify the key aspect that you feel unsure about.

For example, the first step to identify all the description lines in the FASTA file can be decomposed into the sub-steps below.

1. Create a list for storing the description lines
2. Iterate over all the lines in the input file
3. For each line check if it is a description line
4. If it is add the line to the list

Now suppose one felt unsure about how to check if a line was a description line. This would represent a good place to start experimenting with some actual code. Can we come up with a code snippet that can check if a string is a FASTA description line?

## Measuring success¶

Up until this point we have been testing our code snippets by running the code and verifying the results manually. This is terribly inefficient and soon breaks down once you have more than one thing to test.

It is much better to write tests that can be run automatically. Let’s examine this concept. Create a file named fasta_utils.py and add the Python code below to it.

 1 2 3 4 5 6 7 """Module containing utility functions for working with FASTA files.""" def test_is_description_line(): """Test the is_description_line() function.""" print("Testing the is_description_line() function...") test_is_description_line()

The first line is a module level “docstring” and it is used to explain the intent of the module as a whole. A “module” is basically a file with a .py extension, i.e. a Python file. Modules are used to group related functionality together.

Lines three to five represent a function named test_is_description_line(). Python functions consist of two parts: the function definition (line 3) and the function body (lines 4 and 5).

In Python functions are defined using the def keyword. Note that the def keyword is followed by the name of the function. The name of the function is followed by a parenthesized set of arguments, in this case the function takes no arguments. The end of the function definition is marked using a colon.

The body of the function, lines four and five, need to be indented. The standard in Python is to use four white spaces to indent code blocks.

Warning

Whitespace really matters in Python! If your code is not correctly aligned you will see IndentationError messages telling you that everything is not as it should be. You will also run into IndentationError messages if you mix white spaces and tabs.

The first line of the function body, line four, is a docstring explaining the intent of the function. Line five makes use of the built-in print() function to write a string to the standard output stream. Python’s built-in print() function is similar to the echo command we used earlier in Keeping track of your work.

Finally, on line seven the test_is_description_line() function is called, i.e. the logic within the function’s body is executed. In this instance this means that the "Testing the is_description_line() function..." string is written to the standard output stream.

Let us run this script in a terminal. To run a Python script we use the python command followed by the name of the script.

$python fasta_utils.py Testing the is_description_line() function...  So far so good? At the moment our test_is_description_line() function does not actually test anything. Let us rectify that now.  1 2 3 4 5 6 7 8 """Module containing utility functions for working with FASTA files.""" def test_is_description_line(): """Test the is_description_line() function.""" print("Testing the is_description_line() function...") assert is_description_line(">This is a description line") is True test_is_description_line()  There are quite a few things going on in the newly added line. First of all it makes use of three built-in features of Python: the assert and is keywords, as well as the True constant. Let’s work through these in reverse order. Python has some built-in constants, most notably True, False and None. The True and False constants are the only instances of the bool (boolean) type and None is often used to represent the absence of a value. In Python is is an operator that checks for object identity, i.e. if the object returned by the is_description_line() function and True are the same object. If they are the same object the comparison evaluates to True if not it evaluates to False. The assert keyword is used to insert debugging statements into a program. It provides a means to ensure that the state of a program is as expected. If the statement evaluates to False an AssertionError is raised. So, what will happen if we run the code in its current form? Well, we have not yet defined the is_description_line() function, so Python will raise a NameError. Let us run the code. $ python fasta_utils.py
Testing the is_description_line() function...
Traceback (most recent call last):
File "fasta_utils.py", line 8, in <module>
test_is_description_line()
File "fasta_utils.py", line 6, in test_is_description_line
assert is_description_line(">This is a description line") is True
NameError: global name 'is_description_line' is not defined

Great now we are getting somewhere! What? Well, we have impemented some code to test the functionality of the is_description_line() and it tells us that the function does not exist. This is useful information. Let us add a placeholder is_description_line() function to the Python module.

  1 2 3 4 5 6 7 8 9 10 11 """Module containing utility functions for working with FASTA files.""" def is_description_line(line): """Return True if the line is a FASTA description line.""" def test_is_description_line(): """Test the is_description_line() function.""" print("Testing the is_description_line() function...") assert is_description_line(">This is a description line") is True test_is_description_line()

Note that the function we have added on lines three and four currently does nothing. By default the function will return None. However, when we run the script we should no longer get a NameError. Let’s find out what happens when we run the code.

$python fasta_utils.py Testing the is_description_line() function... Traceback (most recent call last): File "fasta_utils.py", line 11, in <module> test_is_description_line() File "fasta_utils.py", line 9, in test_is_description_line assert is_description_line(">This is a description line") is True AssertionError  More progress! Now we see the expected AssertionError, becuase None is not True. Let us add some code to try to get rid of this error message. To achieve this we simply need to make the function return True.   1 2 3 4 5 6 7 8 9 10 11 12 """Module containing utility functions for working with FASTA files.""" def is_description_line(line): """Return True if the line is a FASTA description line.""" return True def test_is_description_line(): """Test the is_description_line() function.""" print("Testing the is_description_line() function...") assert is_description_line(">This is a description line") is True test_is_description_line()  Now, we can run the code again. $ python fasta_utils.py
Testing the is_description_line() function...

No error message, the code is now working to the specification described in the test. However, the test does not specify what the behaviour should be for a biological sequence line. Let us add another assert statement to specify this.

  1 2 3 4 5 6 7 8 9 10 11 12 13 """Module containing utility functions for working with FASTA files.""" def is_description_line(line): """Return True if the line is a FASTA description line.""" return True def test_is_description_line(): """Test the is_description_line() function.""" print("Testing the is_description_line() function...") assert is_description_line(">This is a description line") is True assert is_description_line("ATCG") is False test_is_description_line()

Now we can run the code again.

$python fasta_utils.py Testing the is_description_line() function... Traceback (most recent call last): File "fasta_utils.py", line 13, in <module> test_is_description_line() File "fasta_utils.py", line 11, in test_is_description_line assert is_description_line("ATCG") is False AssertionError  More progress, we now have a test to ensure that the is_description_line() function returns False when the input line is a sequence. Let us try to implement the desired functionality to make the test pass. For this we will use the startswith() method, that is built into strings, to check if the string starts with a greater than (>) character.   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 """Module containing utility functions for working with FASTA files.""" def is_description_line(line): """Return True if the line is a FASTA description line.""" if line.startswith(">"): return True else: return False def test_is_description_line(): """Test the is_description_line() function.""" print("Testing the is_description_line() function...") assert is_description_line(">This is a description line") is True assert is_description_line("ATCG") is False test_is_description_line()  In the code above we make use of conditional logic, i.e. if something is True the code does something otherwise it does something else. As mentioned previously whitespace is important in Python and four spaces are used to indent the lines after the if and else statements to tell Python which statement(s) belong in the conditional code blocks. In this case we only have one statement per conditional, but it is possible to group several statements together based on their indentation. Let us test the code again. $ python fasta_utils.py
Testing the is_description_line() function...

Fantastic the code behaves in the way that we want it to behave!

However, the current implementation of the is_description_line() function is a little bit verbose. Do we really need the else conditional? What would happen if it was not there?

The beauty of tests now become more apparent. We can start experimenting with the implementation of a function and feel confident that we are not breaking it. As long as the tests do not fail that is!

Let us test out our hypothesis that the else conditional is redundant by removing it and de-denting the return False statement.

 1 2 3 4 5 6 7 """Module containing utility functions for working with FASTA files.""" def is_description_line(line): """Return True if the line is a FASTA description line.""" if line.startswith(">"): return True return False

Now we can simply run the tests to see what happens.

$python fasta_utils.py Testing the is_description_line() function...  Amazing, we just made a change to our code and we can feel pretty sure that it is still working as intended. This is very powerful. Incidentally, the code works by virtue of the fact that if the return True statement is reached the function returns True and exits, as such the second return statement is not called. The methodology used in this section is known as Test-Driven Development, often referred to as TDD. It involves three steps: 1. Write a test 2. Write minimal code to make the test pass 3. Refactor the code if necessary In this instance we started off by writing a test checking that the is_description_line() function returned True when the input was a description line. We then added minimal code to make the test pass, i.e. we simply made the function return True. At this point no refactoring was needed so we added another test to check that the function returned False when the input was a sequence line. We then added some naive code to make the tests pass. At this point, we believed that there was scope to improve the implementation of the function, so we refactored it to remove the redundant else statement. Well done! That was a lot of information. Go make yourself a cup of tea. ## More string processing¶ Because information about DNA and proteins are often stored in plain text files many aspects of biological data processing involves manipulating text. In computing text is often referred to as strings of characters. String manipulation is is therefore a common task both for processing biological sequences and for interpreting sequence identifiers. This section will provide a brief summary of how Python can be used for such string processing. Start Python in its interactive mode by typing python into the terminal. $ python
Python 2.7.10 (default, Jul 14 2015, 19:46:27)
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.39)] on darwin
>>>

### The Python string object¶

When parsing in strings from a text file one often has to deal with lines that have leading and/or trailing white spaces. Commonly one wants to get rid of them. This can be achieved using the strip() method built into the string object.

>>> "  text with leading/trailing spaces ".strip()

Another common use case is to replace a word in a line. For example, when we strip out the leading and trailing white spaces one might want to update the word “with” to “without” to make the resulting string reflect its current state. This can be achieved using the replace() method.

>>> "  text with leading/trailing spaces ".strip().replace("with", "without")

Note

In the example above we chain the strip() and replace() methods together. In practise this means that the replace() methods acts on the return value of the strip() method.

Earlier we saw how the startswith() method can be used to identify FASTA description lines.

>>> ">MySeq1|description line".startswith(">")
True

The endswith() method complements the startswith() method and is often used to examine file extensions.

>>> "/home/olsson/images/profile.png".endswith("png")
True

This example above only works if the file extension is in lower case.

>>> "/home/olsson/images/profile.PNG".endswith("png")
False

However, we can overcome this issue by adding a call to the lower() method, which converts the string to lower case.

>>> "/home/olsson/images/profile.PNG".lower().endswith("png")
True

Another common use case is to search for a particular string within another string. For example one might want to find out if the UniProt identifier “Q6GZX4” is present in a FASTA description line. To achieve this one can use the find() method, which returns the index position (zero-based) where the search term was first identified.

>>> ">sp|Q6GZX4|001R_FRG3G".find("Q6GZX4")
4

If the search term is not identified find() returns -1.

>>> ">sp|P31946|1433B_HUMAN".find("Q6GZX4")
-1

When iterating over lines in a file one often wants to split the line based on a delimiter. This can be achieved using the split() method. By default this splits on white space characters and returns a list of strings.

>>> "text without leading/trailing spaces".split()

A different delimiter can be used by providing it as an argument to the split() method.

>>> ">sp|Q6GZX4|001R_FRG3G".split("|")
['>sp', 'Q6GZX4', '001R_FRG3G']

There are many variations on the string operators described above. It is useful to familiarise yourself with the Python documentation on strings.

### Regular expressions¶

Regular expressions can be defined as a series of characters that define a search pattern.

Regular expressions can be very powerful. However, they can be difficult to build up. Often it is a process of trial and error. This means that once they have been created, and the trial and error process has been forgotten, it can be extremely difficult to understand what the regular expression does and why it is constructed the way it is.

Warning

Use regular expression as a last resort. A good rule of thumb is to always try to use string operations to implement the desired functionality and only switch to regular expressions when the code implemented using these operations becomes more difficult to understand than the equivalent regular expression.

To use regular expressions in Python we need to import the re module. The re module is part of Python’s standard library. Importing modules in Python is achieved using the import keyword.

>>> import re

Let us store a FASTA description line in a variable.

>>> fasta_desc = ">sp|Q6GZX4|001R_FRG3G"

Now, let us search for the UniProt identifer Q6GZX4 within the line.

>>> re.search(r"Q6GZX4", fasta_desc)
<_sre.SRE_Match object at 0x...>

There are two things to note here:

1. We use a raw string to represent our regular expression, i.e. the string prefixed with an r
2. The regular expression search() method returns a match object (or None if no match is found)

The index of the first matched character can be accessed using the match object’s start() method. The match object also has an end() method that returns the index of the last character + 1.

>>> match = re.search(r"Q6GZX4", fasta_desc)
>>> if match:
...     print(fasta_desc[match.start():match.end()])
...
Q6GZX4

In the above we make use of the fact that Python strings support slicing. Slicing is a means to access a subsection of a sequence. The [start:end] syntax is inclusive for the start index and exclusive for the end index.

>>> "012345"[2:4]
'23'

To see the merit of regular expressions we need to create one that matches more than one thing. For example a regular expression that could match all the patterns id0, id1, ..., id9.

Now suppose that we had a list containing FASTA description lines with these types of identifiers.

>>> fasta_desc_list = [">id0 match this",
...                    ">id9 and this",
...                    ">id100 but not this (initially)",
...                    "AATCG"]
...

Note that the list above also contains a sequence line that we never want to match.

Let us loop over the items in this list and print out the lines that match our identifier regular expression.

>>> for line in fasta_desc_list:
...     if re.search(r">id[0-9]\s", line):
...         print(line)
...
>id0 match this
>id9 and this

There are two noteworthy aspects of the regular expression. Firstly, the [0-9] syntax means match any digit. Secondly, the \s regular expression meta character means match any white space character.

If one wanted to create a regular expression to match an identifier with an arbitrary number of digits one can make use of the * meta character, which causes the regular expression to match the preceding expression 0 or more times.

>>> for line in fasta_desc_list:
...     if re.search(r">id[0-9]*\s", line):
...         print(line)
...
>id0 match this
>id9 and this
>id100 but not this (initially)

It is possible to extract specific pieces of information from a line using regular expressions. This uses a concept known as “groups”, which are indicated using parenthesis. Let us try to extract the UniProt identifier from a FASTA description line.

>>> print(fasta_desc)
>sp|Q6GZX4|001R_FRG3G
>>> match = re.search(r">sp\|([A-Z,0-9]*)\|", fasta_desc)

Warning

Note how horrible and incomprehensible the regular expression is.

It took me a couple of attempts to get this regular expression right as I forgot that | is a regular expression meta character that needs to be escaped using a backslash \.

The regular expression representing the UniProt idendifier [A-Z,0-9]* means match capital letters (A-Z) and digits (0-9) zero or more times (*). The UniProt regular expression is enclosed in parenthesis. The parenthesis denote that the UniProt identifier is a group that we would like access to. In other words, the purpose of a group is to give the user access to a section of interest within the regular expression.

>>> match.groups()
('Q6GZX4',)
>>> match.group(0)  # Everything matched by the regular expression.
'>sp|Q6GZX4|'
>>> match.group(1)
'Q6GZX4'

Note

There is a difference between the groups() and the group() methods. The former returns a tuple containing all the groups defined in the regular expression. The latter takes an integer as input and returns a specific group. However, confusingly group(0) returns everything matched by the regular expression and group(1) returns the first group making the group() method appear as if it used a one-based indexing scheme.

Finally Let us have a look at a common pitfall when using regular expressions in Python: the difference between the methods search() and match().

>>> print(re.search(r"cat", "my cat has a hat"))
<_sre.SRE_Match object at 0x...>
>>> print(re.match(r"cat", "my cat has a hat"))
None

Basically match() only looks for a match at the beginning of the string to be searched. For more information see the search() vs match() section in the Python documentation.

There is a lot more to regular expressions in particular all the meta characters. For more information have a look at the regular expressions operations section in the Python documentation.

## Extracting the organism name¶

Armed with our new found knowledge of string processing let’s create a function for extracting the organism name from a SwissProt FASTA description line. In other words given the lines:

>sp|P01090|2SS2_BRANA Napin-2 OS=Brassica napus PE=2 SV=2
>sp|Q15942|ZYX_HUMAN Zyxin OS=Homo sapiens GN=ZYX PE=1 SV=1
>sp|Q6QGT3|A1_BPT5 A1 protein OS=Escherichia phage T5 GN=A1 PE=2 SV=1

We would like to extract the strings:

Brassica napus
Homo sapiens
Escherichia phage T5

There are three things which are worth noting:

1. The organism name string is always preceeded by the key OS (Organism Name)
2. The organism name string can contain more than two words
3. The two letter key after the organism name string can vary, in the case above we see both PS (Protein Existence) and GE (Gene Name)

The three FASTA description lines examined above provide an excellent basis for creating a test for the function that we want. Add the lines below to your fasta_utils.py file.

 15 16 17 18 19 20 21 22 23 24 25 def test_extract_organism_name(): """Test the extract_organism_name() function.""" print("Testing the extract_organism_name() function...") lines = [">sp|P01090|2SS2_BRANA Napin-2 OS=Brassica napus PE=2 SV=2", ">sp|Q15942|ZYX_HUMAN Zyxin OS=Homo sapiens GN=ZYX PE=1 SV=1", ">sp|Q6QGT3|A1_BPT5 A1 protein OS=Escherichia phage T5 GN=A1 PE=2 SV=1"] organism_names = ["Brassica napus", "Homo sapiens", "Escherichia phage T5"] for line, organism_name in zip(lines, organism_names): assert extract_organism_name(line) == organism_name test_is_description_line()

In the above we make use of pythons built-in zip() function. This function takes two lists as inputs and returns a list with paired values from the input lists.

>>> zip(["a", "b", "c"], [1, 2, 3])
[('a', 1), ('b', 2), ('c', 3)]

Let’s make sure that the tests fail.

$python fasta_utils.py Testing the is_description_line() function...  What, no error message, what is going on? Ah, we added the test, but forgot to add a line to call it. Let’s rectify that.  25 26 test_is_description_line() test_extract_organism_name()  Let’s try again. $ python fasta_utils.py
Testing the is_description_line() function...
Testing the extract_organism_name() function...
Traceback (most recent call last):
File "fasta_utils.py", line 26, in <module>
test_extract_organism_name()
File "fasta_utils.py", line 23, in test_extract_organism_name
assert extract_organism_name(line) == organism_name
NameError: global name 'extract_organism_name' is not defined

Success! We now have a failing test informing us that we need to create the extract_organism_name() function. Let’s do that.

 15 16 def extract_organism_name(line): """Return the organism name from a FASTA description line."""

Let’s find out where this minimal implementation gets us.

$python fasta_utils.py Testing the is_description_line() function... Testing the extract_organism_name() function... Traceback (most recent call last): File "fasta_utils.py", line 29, in <module> test_extract_organism_name() File "fasta_utils.py", line 26, in test_extract_organism_name assert extract_organism_name(line) == organism_name AssertionError  The test fails as expected. However, since we are looping over many input lines it would be good to get an idea of which test failed. We can achieve this by making use of the fact that we can provide a custom message to be passed to the AssertionError. Let us pass it the input line. Note the addition of the trailing , line in line 26.  25 26  for line, organism_name in zip(lines, organism_names): assert extract_organism_name(line) == organism_name, line  Let’s see what we get now. $ python fasta_utils.py
Testing the is_description_line() function...
Testing the extract_organism_name() function...
Traceback (most recent call last):
File "fasta_utils.py", line 29, in <module>
test_extract_organism_name()
File "fasta_utils.py", line 26, in test_extract_organism_name
assert extract_organism_name(line) == organism_name, line
AssertionError: >sp|P01090|2SS2_BRANA Napin-2 OS=Brassica napus PE=2 SV=2

Much better! Let us try to implement a basic regular expression to make this pass. First of all we need to make sure we import the re module.

 1 2 3 """Module containing utility functions for working with FASTA files.""" import re

Then we can implement a regular expression to try to extract the organism name.

 17 18 19 20 def extract_organism_name(line): """Return the organism name from a FASTA description line.""" match = re.search(r"OS=(.*) PE=", line) return match.group(1)

Remember the .* means match any character zero or more times and the surrounding parenthesis creates a group.

Let us see what happens now.

$python fasta_utils.py Testing the is_description_line() function... Testing the extract_organism_name() function... Traceback (most recent call last): File "fasta_utils.py", line 34, in <module> test_extract_organism_name() File "fasta_utils.py", line 31, in test_extract_organism_name assert extract_organism_name(line) == organism_name, line AssertionError: >sp|Q15942|ZYX_HUMAN Zyxin OS=Homo sapiens GN=ZYX PE=1 SV=1  Progress! We are now seeing a different error message. The issue is that the key after the regular expression is GN rather than PE. Let us try to rectify that.  17 18 19 20 def extract_organism_name(line): """Return the organsim name from a FASTA description line.""" match = re.search(r"OS=(.*) [A-Z]{2}=", line) return match.group(1)  The regular expression now states that instead of PE it wants any capital letter [A-Z] repeated twice {2}. Let’s find out if this fixes the issue. $ python fasta_utils.py
Testing the is_description_line() function...
Testing the extract_organism_name() function...
Traceback (most recent call last):
File "fasta_utils.py", line 33, in <module>
test_extract_organism_name()
File "fasta_utils.py", line 30, in test_extract_organism_name
assert extract_organism_name(line) == organism_name, line
AssertionError: >sp|P01090|2SS2_BRANA Napin-2 OS=Brassica napus PE=2 SV=2

What, back at square one again? As mentioned previously, regular expressions can be painful and should only be used as a last resort. This also exemplifies why it is important to have tests. Sometimes you think you make an innocuous change, but instead things just fall apart.

At this stage the error message is not very useful, let us change it to print out the value returned by the function instead.

 29 30  for line, organism_name in zip(lines, organism_names): assert extract_organism_name(line) == organism_name, extract_organism_name(line)

Now, let’s see what is going on.

$python fasta_utils.py Testing the is_description_line() function... Testing the extract_organism_name() function... Traceback (most recent call last): File "fasta_utils.py", line 33, in <module> test_extract_organism_name() File "fasta_utils.py", line 30, in test_extract_organism_name assert extract_organism_name(line) == organism_name, extract_organism_name(line) AssertionError: Brassica napus PE=2  Our regular expression is basically matching too much. The reason for this is that the * meta character acts in a “greedy” fashion matching as much as possible, see Fig. 20. In this case the PE=2 is included in the match group as [A-Z]{2} is matched by the SV= key at the end of the line. The fix is to make the * meta character act in a “non-greedy” fashion. This is achieved by adding a ? suffix to it.  17 18 19 20 def extract_organism_name(line): """Return the organism name from a FASTA description line.""" match = re.search(r"OS=(.*?) [A-Z]{2}=", line) return match.group(1)  Let’s find out what happens now. $ python fasta_utils.py
Testing the is_description_line() function...
Testing the extract_organism_name() function...

All the tests pass! Well done, time for another cup of tea.

## Only running tests when the module is called directly¶

In Python, modules provide a means to group related functionality together. For example we have already looked at and made use of the re module, which groups functionality for working with regular expressions.

In Python any file with a .py extension is a module. This means that the file that we have been creating, fasta_utils.py, is a module.

To make use of the functionality within a module one needs to import it. Let’s try this out in an interactive session

$python  Now we can import the module. >>> import fasta_utils Testing the is_description_line() function... Testing the extract_organism_name() function...  Note that the tests run just like when we call the fasta_utils.py script directly. This is an undesired side effect of the current implementation. It would be better if the tests were not run when the module is imported. To improve the behaviour of the fasta_utils module we will make use of a special Python attribute called __name__, which provides a string representation of scope. When commands are run from a script or the interactive prompt the name attribute is set to __main__. When a module is imported the __name__ attribute is set to the name of the module. >>> print(__name__) __main__ >>> print(fasta_utils.__name__) fasta_utils  Using this information we can update the fasta_utils.py file with the changes highlighted below.  32 33 34 if __name__ == "__main__": test_is_description_line() test_extract_organism_name()  Let us make sure that the tests still run if we run the script directly. $ python fasta_utils.py
Testing the is_description_line() function...
Testing the extract_organism_name() function...

Now we can reload the module in the interactive prompt we were working in earlier to make sure that the tests no longer get executed.

>>> reload(fasta_utils)
<module 'fasta_utils' from 'fasta_utils.py'>

The test suite is no longer being called upon loading the module and no lines logging the progress of the testing suite is being printed to the terminal.

Note that simply calling the import fasta_utils command again will not actually detect the changes that we made to the fasta_utils.py file, which is why we make use of Python’s built-in reload() function. Alternatively, one could have exited the Python shell, using Ctrl-D or the exit() function, and then started a new interactive Python session and imported the fasta_utils module again.

## Counting the number of unique organisms¶

We can now use the fasta_utils module to start answering some of the biological questions that we posed at the beginning of this chapter. For now let us do this using an interactive Python shell.

$python Python 2.7.10 (default, Jul 14 2015, 19:46:27) [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.39)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>>  Now we will start by importing the modules that we want to use. In this case fasta_utils for processing the data and gzip for opening and reading in the data from the uniprot_sprot.2015-11-26.fasta.gz file. >>> import fasta_utils >>> import gzip  Now we create a list for storing all the FASTA description lines. >>> fasta_desc_lines = []  We can then open the file using the gzip.open() function. Note that this returns a file handle. >>> file_handle = gzip.open("../data/uniprot_sprot.2015-11-26.fasta.gz")  Using a for loop we can iterate over all the lines in the input file. >>> for line in file_handle: ... if fasta_utils.is_description_line(line): ... fasta_desc_lines.append(line) ...  When we are finished with the input file we must remember to close it. >>> file_handle.close()  Let’s check that the number of FASTA description lines using the built-in len() function. This function returns the number of items in the list, i.e. the length of the list. >>> len(fasta_desc_lines) 549832  Okay now it is time to find the number of unique organisms. For this we will make use of a data structure called set. In Python sets are used to compare collections of unique elements. This means that sets are ideally suited for operations that you may associate with Venn diagrams. However, in this instance we simply use the set data structure to ensure that we only get one unique representative of each organism. In other words even if one calls the set.add() function several times with the same item the item will only occur once in the set. >>> organisms = set() >>> for line in fasta_desc_lines: ... s = fasta_utils.extract_organism_name(line) ... organisms.add(s) ... >>> len(organisms) 13251  Great, now we know that there are 13,251 unique organisms represented in the FASTA file. ## Finding the number of variants per species and the number of proteins per variant¶ Suppose we had a FASTA file containing only four entries with the description lines below. >sp|P12334|AZUR1_METJ Azurin iso-1 OS=Methylomonas sp. (strain J) PE=1 SV=2 >sp|P12335|AZUR2_METJ Azurin iso-2 OS=Methylomonas sp. (strain J) PE=1 SV=1 >sp|P23827|ECOT_ECOLI Ecotin OS=Escherichia coli (strain K12) GN=eco PE=1 SV=1 >sp|B6I1A7|ECOT_ECOSE Ecotin OS=Escherichia coli (strain SE11) GN=eco PE=3 SV=1  The analysis in the previous section would have identified these as three separate entities. Methylomonas sp. (strain J) Escherichia coli (strain K12) Escherichia coli (strain SE11)  Now, suppose that we wanted to find out how many variants there were of each species. In the example above there would be be one variant of Methylomonas sp. and two variants of Escherichia coli. Furthermore, suppose that we also wanted to find out how many proteins were associated with each variant. We could achieve this by creating a nested data structure using Python’s built in dictionary type. At the top level we should have a dictionary whose keys were the species, e.g. Escherichia coli. The values in the top level dictionary should themselves be dictionaries. The keys of the nested dictionaries should be the full organism name, e.g. Escherichia coli (strain K12). The values in the nested dictionary should be an integer representing the number of proteins found for that organism. Below is a YAML representation of the data structure that should be created from the four entries above. --- Methylomonas sp.: Methylomonas sp. (strain J): 2 Escherichia coli: Escherichia coli (strain K12): 1 Escherichia coli (strain SE11): 1  So what type of functionality would we need to achieve this? First of all we need a function that given an organism name returns the associated species. In other words something that converts Escherichia coli (strain K12) to Escherichia coli. Secondly, we need a function that given a list of organism names returns the data structure described above. Let us start by creating a test for converting the organism name into a species. Add the test below to fasta_utils.py.  32 33 34 35 36 37 38 39 40 def test_organism_name2species(): print("Testing the organism_name2species() function...") assert organism_name2species("Methylomonas sp. (strain J)") == "Methylomonas sp." assert organism_name2species("Homo sapiens") == "Homo sapiens" if __name__ == "__main__": test_is_description_line() test_extract_organism_name() test_organism_name2species()  Let’s find out what error this gives us. $ python fasta_utils.py
Testing the is_description_line() function...
Testing the extract_organism_name() function...
Testing the organism_name2species() function...
Traceback (most recent call last):
File "fasta_utils.py", line 40, in <module>
test_organism_name2species()
File "fasta_utils.py", line 34, in test_organism_name2species
assert organism_name2species("Methylomonas sp. (strain J)") == "Methylomonas sp."
NameError: global name 'organism_name2species' is not defined

The error message is telling us that we need to define the organism_name2species() function. Add the lines below to define it.

 32 33 def organism_name2species(organism_name): """Return the species from the FASTA organism name."""

Now we get a new error message when we run the tests.

$python fasta_utils.py Testing the is_description_line() function... Testing the extract_organism_name() function... Testing the organism_name2species() function... Traceback (most recent call last): File "fasta_utils.py", line 43, in <module> test_organism_name2species() File "fasta_utils.py", line 37, in test_organism_name2species assert organism_name2species("Methylomonas sp. (strain J)") == "Methylomonas sp." AssertionError  Great, let us add some logic to the function.  32 33 34 35 def organism_name2species(organism_name): """Return the species from the FASTA organism name.""" words = organism_name.split() return words[0] + " " + words[1]  Above, we split the organism name based on whitespace separators and return the first two words joined by a space character. Note In Python, and many other scripting languages, strings can be concatenated using the + operator. >>> "hello" + " " + "world" 'hello world'  Time to test the code again. $ python fasta_utils.py
Testing the is_description_line() function...
Testing the extract_organism_name() function...
Testing the organism_name2species() function...

Great, the function is working! Let us define a new test for the function that will generate the nested data structure we described earlier.

 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 def test_summarise_species_protein_data(): print("Testing summarise_species_protein_data() function...") fasta_desc_lines = [ ">sp|P12334|AZUR1_METJ Azurin iso-1 OS=Methylomonas sp. (strain J) PE=1 SV=2", ">sp|P12335|AZUR2_METJ Azurin iso-2 OS=Methylomonas sp. (strain J) PE=1 SV=1", ">sp|P23827|ECOT_ECOLI Ecotin OS=Escherichia coli (strain K12) GN=eco PE=1 SV=1", ">sp|B6I1A7|ECOT_ECOSE Ecotin OS=Escherichia coli (strain SE11) GN=eco PE=3 SV=1" ] summary = summarise_species_protein_data(fasta_desc_lines) # The top level dictionary will contain two entries. assert len(summary) == 2 assert "Methylomonas sp." in summary assert "Escherichia coli" in summary # The value of the Methylomonas sp. entry is a dictionary with one # entry in it. assert len(summary["Methylomonas sp."]) == 1 assert summary["Methylomonas sp."]["Methylomonas sp. (strain J)"] == 2 # The value of the Escherichia coli entry is a dictionary with two # entries in it. assert len(summary["Escherichia coli"]) == 2 assert summary["Escherichia coli"]["Escherichia coli (strain K12)"] == 1 assert summary["Escherichia coli"]["Escherichia coli (strain SE11)"] == 1 if __name__ == "__main__": test_is_description_line() test_extract_organism_name() test_organism_name2species() test_summarise_species_protein_data()

This should all be getting familiar now. Time to run the tests.

$python fasta_utils.py Testing the is_description_line() function... Testing the extract_organism_name() function... Testing the organism_name2species() function... Testing summarise_species_protein_data() function... Traceback (most recent call last): File "fasta_utils.py", line 70, in <module> test_summarise_species_protein_data() File "fasta_utils.py", line 50, in test_summarise_species_protein_data summary = summarise_species_protein_data(fasta_desc_lines) NameError: global name 'summarise_species_protein_data' is not defined  Again we start by defining the function.  42 43 def summarise_species_protein_data(fasta_desc_lines): """Return data structure summarising the organism and protein data."""  And then we run the tests again. $ python fasta_utils.py
Testing the is_description_line() function...
Testing the extract_organism_name() function...
Testing the organism_name2species() function...
Testing summarise_species_protein_data() function...
Traceback (most recent call last):
File "fasta_utils.py", line 74, in <module>
test_summarise_species_protein_data()
File "fasta_utils.py", line 57, in test_summarise_species_protein_data
assert len(summary) == 2
TypeError: object of type 'NoneType' has no len()

 42 43 44 45 46 47 48 49 50 51 def summarise_species_protein_data(fasta_desc_lines): """Return data structure summarising the organism and protein data.""" summary = dict() for line in fasta_desc_lines: variant_name = extract_organism_name(line) species_name = organism_name2species(variant_name) variant_dict = summary.get(species_name, dict()) variant_dict[variant_name] = variant_dict.get(variant_name, 0) + 1 summary[species_name] = variant_dict return summary

In the above we make use of the dictionary’s built-in get() method, which returns the value associated with the key provided as the first argument. If the key does not yet exist in the dictionary it returns the second argument, the default value.

On line 48 we try to access a variant_dict dictionary from within the summary dictionary, we are in other words working with a nested data structure. If there is no entry associated with the species_name the get() method will return an empty dictionary (dict()).

On line 49 we keep count of the number of times a variant of a species, as defined by the full organism name, has been observed. In this instance we use the get() method to get the number of times the variant has been observed previously and then we add 1 to it. If the variant has never been observed previously it will not be in the dictionary and the variant_dict.get(variant_name, 0) method call will return 0 (the default value specified).

Line 50 updates the top level dictionary by adding the variant_dict dictionary to the summary dictionary.

Let’s see if it the implementation works as expected.

$python fasta_utils.py Testing the is_description_line() function... Testing the extract_organism_name() function... Testing the organism_name2species() function... Testing summarise_species_protein_data() function...  Hurray! Finally, let us write a separate script to convert an input FASTA file into a YAML summary file. Create the file fasta2yaml_summary.py and add the code below to it. #!/usr/bin/env python import sys import yaml import fasta_utils fasta_desc_lines = list() with sys.stdin as fh: for line in fh: if fasta_utils.is_description_line(line): fasta_desc_lines.append(line) summary = fasta_utils.summarise_species_protein_data(fasta_desc_lines) yaml_text = yaml.dump(summary, explicit_start=True, default_flow_style=False) with sys.stdout as fh: fh.write(yaml_text)  In the code above we make use of the yaml module to convert our data structure to the YAML file format. The PyYAML package is not part of the Python’s standard library, but it is easily installed using pip (for more detailed instructions see Installing Python packages). $ sudo pip install pyyaml

The script also makes use of sys.stdin and sys.stdout to read from the standard input stream and write to the standard output stream, respectively. This means that we can pipe in the content to our script and pipe output from our script. For example to examine the YAML output using the less pager one could use the command below.

\$ gunzip -c data/uniprot_sprot.2015-11-26.fasta.gz | python fasta2yaml_summary.py | less

This immediately reveals that there are organisms in the SwissProt FASTA file that have few protein associated with them.

Note

Remember the benefits of using pipes was described in Creating a work flow using pipes.

---
AKT8 murine:
AKT8 murine leukemia virus: 1
AKV murine:
AKV murine leukemia virus: 3
Abelson murine:
Abelson murine leukemia virus: 3
Abies alba:
Abies alba: 1
Abies balsamea:
Abies balsamea: 3

Note

The formatting in the YAML file above was created for us by the call to the yaml.dump() method in the fasta2yaml_summary.py script.

Great work!

## Key concepts¶

• Think about the problem at hand before you start coding
• Python has many built-in packages for performing complex tasks, in this chapter we used the re package for working with regular expressions
• There are also many third party Python packages that can be installed, in this chapter we made use of the yaml package for writing out a data structure as a YAML file


text:
## Chapter 3: Exercises

1. (*)About k-means clustering: Given a set of $n$ points in a $d$-dim space $\mathbf{X}=\{\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_n\}$, we want to find another set of $m$ points in the same space $\mathbf{C}=\{\mathbf{c}_1, \mathbf{c}_2, \dots, \mathbf{c}_m\}$, with $m \lt \lt n$.
1. What is the goal of k-means clustering?
2. Define the objective function of k-means clustering. (Please explain your math notations in detail.)
3. Describe how the objective function can be minimized in a step-by-step manner.
4. Why the objective function is guaranteed to be non-increasing during the iterations?
2. (*)K-means clustering for image data compression: Please explain briefly how k-means clustering can be used for image data compression.
3. (*)Some identities for k-means clustering: Here are some identities for k-means and k-medians clustering.
1. Given a vector $\mathbf{x}=[x_1, x_2, \dots, x_n]$, prove that $\arg \min_s \sum_{i=1}^n (s-x_i)^2 = mean(\mathbf{x})$.
2. Given a vector $\mathbf{x}=[x_1, x_2, \dots, x_n]$, prove that $\arg \min_s \sum_{i=1}^n |s-x_i| = median(\mathbf{x})$.
3. Given a set of vectors $\mathbf{X}=[\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_n]$, prove that $\arg \min_\mathbf{s} \sum_{i=1}^n \|\mathbf{s}-\mathbf{x}_i\|_2^2 = mean(\mathbf{X}, 2)$.
4. Given a set of vectors $\mathbf{X}=[\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_n]$, prove that $\arg \min_\mathbf{s} \sum_{i=1}^n \|\mathbf{s}-\mathbf{x}_i\|_1 = median(\mathbf{X}, 2)$.
Note that in the above expressions, $\|\mathbf{x}\|_p=\left(\sum_{i=1}^n x_i^p\right)^{1/p}$.
4. (*) Objective function w.r.t. the number of clusters: Write a script that obtains the dataset via the following statement: DS = dcData(5); Then use kMeansClustering.m for k-means clustering, according to the following criterion.
1. Let the number of clusters start from 1 to 10 and plot the objective function. The plot should be similar to the one shown next:
2. From the scatter plot of the data set (see this example), we know the best number of clusters is 4 by visual inspection. Can you think of a reliable way to determine the best number of clusters automatically? Please implement your method in the script and print out the best number of clusters directly. This is your first attemp at cluster validation which uses a certain criterion to determine the best number of clusters automatically.
3. To see how your method for cluster validation works, try the following datasets:
• DS = dcData(1);
• DS = dcData(3);
• DS = dcData(4);
• DS = dcData(6);
Since the datasets are generated by random number, you should put your method in a for loop and run it 10 times to get the average result.
(Hint: try "dcData" to view the scatter plots directly.)
5. (**) Selection of the initial centers: Please modify the first two lines of this example as follows: DS = dcData(1); centerNum=6; After running the script several times, you will notice that the initial centers play an important roles in determining the final clustering result. Since kMeansClustering.m can also take a set of initial centers for clustering, your task is to write an m-file function myVqInitCenter.m that can return a set of centers for better clustering for k-means, with the same usage as vqInitCenter.m. As described in the text, some of the methods include:
1. Randomly select k data points from the dataset.
2. Select the farthest k data points from the mean of the dataset.
3. Select the nearest k data points from the mean of the dataset.
4. Select k data points that have the largest sum of pairwise square distance.
Please implement the last method to see if you can get better results of k-means clustering. In particular, the search procedure is like this:
1. Find two points that has the maximal square distance among all data set. Put these two data set in the pool of centers.
2. Add a new points to the pool that can increase the sum of square distance most.
3. Repeat the previous step until we have k centers.
Use the following statements to test your program: centerNum=6; runNum=20; dist1=zeros(runNum,1); dist2=zeros(runNum,1); % ====== Get initial cluster centers for i=1:runNum DS = dcData(1); initCenter=vqInitCenter(DS.input, centerNum, 1); [center, U, distortion, allCenters] = kMeansClustering(DS.input, initCenter); dist1(i)=distortion(end); initCenter=myVqInitCenter(DS.input, centerNum); [center, U, distortion, allCenters] = kMeansClustering(DS.input, initCenter); dist2(i)=distortion(end); end plot(1:runNum, dist1, 'o-', 1:runNum, dist2, 'o-'); xlabel('Run index'); ylabel('Distortion'); legend('vqInitCenter', 'myVqInitCenter'); Your plot should be similar to the one shown next:
Hint: You can use the following functions if necessary:
• distSqrPairwise.m (in Machine Learning Toolbox): Compute the pairwise square distance between any two points in a dataset.
• maxxy.m (in Utility Toolbox): Find the maximal element in a 2D matrix.
6. (*) Image compression ratio achieved by k-means clustering: If we use 256 colors to represent an indexed-color image of m by n, what is the compression ratios of the following cases when compared with a true-color image of the same size?
1. Use pixels for k-means clustering.
2. Use 2x2 blocks for k-means clustering.
3. Use 4x4 blocks for k-means clustering.
4. Use qxq blocks for k-means clustering, assuming q is a factor of both m and n.
7. (***) Use k-means clustering on pixels of images for image compression: A true-color image of size mxn is represented by mxn pixels, each consists of 24 bits of RGB colors. On the other hand, each pixel of an mxn index-color image is represented by an p-bit unsigned integer, which serves as an index into a color map of size 2^p by 3. In order to save storage space, we can use k-means clustering to convert a true-color image into an index-color image. The original true-color image requires 24*m*n bits of storage. After the conversion, the index-color image requires p*m*n+8*2^p*3 bits of storage. If we use 256 colors to represent an indexed image of 480*640, the compression ratio is (24*480*640)/(8*480*640+8*256*3) = 2.9925.

The goal of this exercise is to display the images after data compression using k-means clustering. The original true-color image is shown next:

Your mission is to convert the original true-color image into an index-color one using k-means clustering on each pixel represented as a vector of 3 elements of RGB intensity. Let k take the values of 2, 4, 8, 16, 32, 64 and display the corresponding index-color images as shown next. How long does it take to run each k-means clustering? What are the distortion of each images?

Hint:
• You need to read the original image and convert it into a 3x307200 data matrix of "double", where each column is the RGB vector of an pixel.
• Use imread to read an image, image to display an image, and colorbar to display the color map.
• Use kMeansClustering for k-means clustering.
• Use distPairwise to find the distance between two sets of vectors.
8. (***) Use k-means clustering on blocks of images for image compression: Repeat the previous exercise but use blocks of different sizes for k-means clustering.
1. 2x2 blocks (4 pixels, so the dimension is 3*4=12 for k-means clustering). My results are as follows:

2. 4x4 blocks (16 pixels, so the dimension is 3*16=48 for k-means clustering). My results are as follows:

3. 8x8 blocks (64 pixels, so the dimension is 3*64=192 for k-means clustering). My results are as follows:

1. What are the compression rate formula for using blocks of qxq?
2. Which block size give the best compression rate with similar image quality of pixel-based k-means clustering with no. of clusters equal to 64 (image shown next)? Can you show them side-by-side?

9. (***) Use k-means clustering on blocks of images for image compression, with 3 codebooks for R, G, and B: Repeat the previous exercise but use 3 codebooks for R, G, and B. In other words, each layer (R, G, or B) has its own codebook.
10. (***) Image compression function with k-means clustering: Write a function myImageCompress.m that can perform impage compression on blocks of 8x8, with given number of centers/clusters for k-means clustering. The ussage of the function is as follows:
outImage=myImageCompress(inImage, k)
, where
• inImage: input true-color image, with dimension mxnx3.
• k: number of clusters (or clusters)
• outImage: output image of mxnx3, where each 8x8 block has already been replaced by the corresponding cluster center obtained via k-means clustering on 8x8 blocks.
Here are two test cases on the test image:

• Case 1:

Note that I've tried to put as comments as possible in the example. If you still have questions, please contact TA as soon as you can.)

Example 1: myImageCompress01.m% You need to change the following two lines to indicate where you put the Utility Toolbox and Machine Learning Toolbox addpath /users/jang/matlab/toolbox/utility addpath /users/jang/matlab/toolbox/machineLearning close all; % Close all windows clear all; % Cear all variables in the workspace inImage = imread('vivian_orig.png'); % Read image file [m, n]=size(inImage); % For the size of the input image blockSide=8; % For computing compression ratio centerCount=64; % No. of centers/clusters outImage=myImageCompress(inImage, centerCount); % Perform image compression using k-means clustering absDiff=abs(double(inImage)-outImage); % You need to use double() to convert the image data back to double precision distortion=sum(absDiff(:).^2); % Total distortion fprintf('distortion = %f\n', distortion); % Print total distortion fprintf('Average deviation = %f\n', mean(absDiff(:))); % Print average deviation for each pixel component %% Plot the images figure; subplot(211); image(inImage); axis image; title('Original image'); subplot(212); image(uint8(outImage)); axis image; % Display the output image. You need to use uint8() to convert double precision back to unsigned 8-bit integers compressionRatio=8*3*m*n/(log2(centerCount)*(m/blockSide)*(n/blockSide)+8*centerCount*3*blockSide*blockSide); % Compute the compression ratio titleStr=sprintf('Codebook size=%g, block=%dx%d, compression ratio=%g', centerCount, blockSide, blockSide, compressionRatio); % Note that codebook size is equal to the number of centers/clusters title(titleStr);distortion = 69168242.400052 Average deviation = 5.244644

• Case 2:

Example 2: myImageCompress02.maddpath /users/jang/matlab/toolbox/utility addpath /users/jang/matlab/toolbox/machineLearning close all; clear all; inImage = imread('vivian_orig.png'); % Read image [m, n]=size(inImage); % For computing compression ratio blockSide=8; % For computing compression ratio maxI=8; % The max no. of centers is 2^maxI for i=1:maxI tic centerCount=2^i; % No. of centers fprintf('i=%d/%d: no. of centers=%d\n', i, maxI, centerCount); outImage=myImageCompress(inImage, centerCount); %% You can un-comment the following lines to display the resultant images at each iteration % figure; image(uint8(outImage)); axis image; % Display the output image % compressionRatio=8*3*m*n/(log2(centerCount)*(m/blockSide)*(n/blockSide)+8*centerCount*3*blockSide*blockSide); % titleStr=sprintf('Codebook size=%g, block=%dx%d, compression ratio=%g', centerCount, blockSide, blockSide, compressionRatio); % title(titleStr); outputFile=sprintf('vivian_%.3d_%dx%d.png', centerCount, blockSide, blockSide); imwrite(uint8(outImage), outputFile); % Write output image file fprintf('\tTime=%f sec\n', toc); end fileList=dir('vivian*.png'); montage({fileList.name}); % Display all images in a single ploti=1/8: no. of centers=2 Time=0.571409 sec i=2/8: no. of centers=4 Time=0.462520 sec i=3/8: no. of centers=8 Time=0.566101 sec i=4/8: no. of centers=16 Time=0.964025 sec i=5/8: no. of centers=32 Time=1.089227 sec i=6/8: no. of centers=64 Time=1.876587 sec i=7/8: no. of centers=128 Time=2.126622 sec i=8/8: no. of centers=256 Time=3.390745 sec [Warning: Image is too big to fit on screen; displaying at 67%] [> In images.internal.initSize (line 71) In imshow (line 332) In montage (line 156) In myImageCompress02 (line 24) In dummyFunction', 'd:\users\jang\books\goWriteOutputFile.m', 85)" style="font-weight:bold">goWriteOutputFile>dummyFunction (line 85) In goWriteOutputFile (line 55)]

Data Clustering and Pattern Recognition (¸ê®Æ¤À¸s»P¼Ë¦¡¿ë»{)


text:
### Home > PC3 > Chapter 8 > Lesson 8.2.2 > Problem8-77

8-77.
1. Use the properties of logarithms to write each of the following expressions using a single logarithm. Homework Help ✎

1. 2 log(m) – 3 log(n) + log(p)

$\log(m^2)- \log(n^3) + \log(p^{1/2})$

Work inside the parentheses first to make a single logarithm. Then apply the 1/2 as an exponent.


text:
Is there an explicit formula for $\sum x^n/(n!)^2$? - MathOverflow most recent 30 from http://mathoverflow.net 2013-06-19T15:13:21Z http://mathoverflow.net/feeds/question/94119 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/94119/is-there-an-explicit-formula-for-sum-xn-n2 Is there an explicit formula for $\sum x^n/(n!)^2$? xmchenhit 2012-04-15T14:44:41Z 2012-04-16T00:27:23Z
Let $f(x)=\sum_{n\geq 0}\frac{1}{n!n!}x^{n}$, is there an explicit formula for $f(x)$?
 http://mathoverflow.net/questions/94119/is-there-an-explicit-formula-for-sum-xn-n2/94123#94123 Answer by Woett for Is there an explicit formula for $\sum x^n/(n!)^2$? Woett 2012-04-15T15:09:47Z 2012-04-15T16:53:38Z
This is, for $x \ge 0$, a special case of the modified Bessel function of the first kind. Have a look here: http://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html $$\sum_{n = 0}^{\infty} \frac{x^{n}}{(n!)^{2}} = \mathrm{I}_0 \left(2 \sqrt{x}\right)$$


text:
# [NTG-context] Generalized \installcompoundcharacter

Matija Pretnar ntg-context@ntg.nl
Thu, 15 Apr 2004 16:17:02 +0200

With the advent of UTF-8, there comes a great possibility to improve
For example: instead of \int_a^b, you could write S_a^b, where S is the
character for integral. Or greek letters and any other math symbol
imaginable.
I tried using \installcompoundcharacter S {\int} and similar commands
and all the characters show up nicely, but the spacing is all wrong,
because they are interpreted as variables instead of operations,
relations, ...
I think this replacement would be easy to implement (if it isn't
already implemented), because with math symbols, you don't have to
worry about hyphenation. You could simply replace the characters with


text:
# INFRARED EMISSION AS STUDIED THROUGH THE TEMPERATURE DEPENDENCE OF INTEGRATED INTENSITIES OF A MOLECULAR CRYSTAL.

Please use this identifier to cite or link to this item: http://hdl.handle.net/1811/15290

Files Size Format View
1967-N-04.jpg 115.9Kb JPEG image

 Title: INFRARED EMISSION AS STUDIED THROUGH THE TEMPERATURE DEPENDENCE OF INTEGRATED INTENSITIES OF A MOLECULAR CRYSTAL. Creators: Leles, Byron K.; Neff, Vernon D. Issue Date: 1967 Publisher: Ohio State University Abstract: The expression for the integrated intensity of an infrared absorption band is normally obtained by neglecting the probability of spontaneous emission. At high temperatures, this factor cannot be neglected, and in fact gives rise to a decrease in intensity with increasing temperature. We have studied this effect experimentally with the molecular crystal, copper phthalocyanine at temperatures up to $600^{\circ}C$. Expressions based on the assumption of a harmonic potential function have been derived for the temperature dependence of the fundamentals and certain overtones and combinations. These expressions compare favorably with the experimental data. Description: Author Institution: Department of Chemistry, Kent State University URI: http://hdl.handle.net/1811/15290 Other Identifiers: 1967-N-4


text:
GCSE Maths Geometry and Measure

2D Shapes

CSA

Arc Of A Circle

# Arc Of A Circle

Here we will learn about the arc of a circle, including how to find the arc length and work backwards given the arc length to solve problems.

There are also arc of a circle worksheets based on Edexcel, AQA and OCR exam questions, along with further guidance on where to go next if you’re still stuck.

## What is an arc of a circle?

An arc of a circle is a portion of the circumference of a circle.

In GCSE mathematics the arc of a circle can be used to define angles and lengths.

The arc of a circle would be the curved edge of a sector or segment.

There are circle theorems that are defined using arcs of the circle.

For example, \theta is the angle subtended by arc AC.

We can find the length of an arc if we know the angle of the sector.

### Calculating the arc length

The arc length is a fraction of the circumference of the circle.

The fraction will be \frac{\theta}{360}, if \theta is the angle of the sector.

The arc length formula is,

\text{Arc length } = \frac{\theta}{360} \times \pi\times d

\theta = \text{ angle of the sector}

d = \text{ diameter of the circle}

Or

\text{Arc length } = \frac{\theta}{360} \times 2\times\pi \times r

\theta = \text{ angle of the sector}

r = \text{ radius of the circle}

## How to calculate the arc length

In order to calculate the arc length:

1. Find the length of the radius.
2. Find the size of the angle creating the sector.
3. Substitute the values you know into the formula for the arc length.

## Calculating the arc length examples

### Example 1: calculate the arc length given the angle

Calculate the arc length of the sector shown below.

1. Find the length of the radius/diameter.

2Find the size of the angle creating the sector.

Angle = 115^{\circ}

3Substitute the values you know into the formula for the arc length.

As you know the radius you can use the formula which has ‘r’ as a variable.

\begin{aligned} \text{Arc length } &=\frac{\theta}{360} \times 2\times\pi \times r \\\\ &=\frac{115}{360} \times 2\times\pi \times 8 \\\\ &=\frac{46}{9} \pi \end{aligned}

\begin{aligned} \text{Arc length } &=\frac{46}{9}\pi \ cm \\\\ &=16.0570...cm \end{aligned}

The arc measures =16.1 \ cm .

### Example 2: calculate the arc length not given the angle

Calculate the arc length of the sector shown below.

Although we can use the radius or diameter as long as we use the correct formula for each.

The angle creating the sector is the angle included between the two radii.

To find it in this case we need to use the angle fact, angles around a point add to 360^{\circ}.

360-240=120

The angle creating the sector is 120^{\circ}.

\begin{aligned} \text { Arc length }&=\frac{\theta}{360} \times 2 \times \pi \times r \\\\ &=\frac{120}{360} \times 2 \times \pi \times 4 \end{aligned}

\begin{aligned} &=8.37758041 \\\\ &=8.38 \mathrm{~cm} \end{aligned}

### Example 3: finding the arc length not given the radius or angle

Below is a sector of a circle.

The diameter of the circle is 12 \ cm.

Calculate the arc length.

Diameter = 12 \ cm

The radius is half of the diameter, so the radius = 6 \ cm.

Although we can use the radius or diameter as long as we use the correct formula for each.

The angle creating the sector is the angle included between the two radii.

To find it in this case we need to use the angle fact, angles around a point add to 360^{\circ}.

360-205=155

The angle creating the sector is 155^{\circ}.

\begin{aligned} \text { Arc length }&=\frac{\theta}{360} \times 2 \times \pi \times r \\\\ &=\frac{155}{360} \times 2 \times \pi \times 6 \end{aligned}

\begin{aligned} &=16.23156204 \\\\ &=16.2 \mathrm{~cm} \end{aligned}

### Common misconceptions

• Finding the length of the circumference not the length of the arc

Remember to find the fraction of the circle that makes the arc not just the circumference of the whole circle.

For example,

Use \frac{\theta}{360} \times \pi\times d

or \frac{\theta}{360} \times 2\times\pi \times r .

Not \pi\times d

or 2\times\pi \times r.

• Incorrect angle unit set on calculator

At GCSE all angles are measured in degrees. Make sure that your calculator has a small ‘d’ for degrees at the top of the screen rather than an ‘r’ for radians – these are not used until A Level.

### Practice arc of a circle questions

1. What part of a circle does the red line show?

Segment

Major arc

Minor arc

Sector

Part of the circumference is called an arc. The red line is more than half of the circumference, so it is a major arc.

2. The diagram shows a semicircle. What is the length of the curved edge, in terms of \pi.

8 \pi \ cm

8\pi+8 \ cm

4\pi \ cm

2 \pi \ cm

The length of the curved edge of a semicircle can be found by calculating \pi d \div 2.

3. Find the length of the arc of the sector. Give your answer to 1 decimal place.

9.4 \ cm

28.3 \ cm

18.8 \ cm

37.7 \ cm

The length of the curved edge of a quarter-circle can be found by calculating \pi d \div 4.

4. Find the length of the arc of the sector. Give your answer to 1 decimal place.

2.4 \ cm

7.3 \ cm

17.1 \ cm

4.9 \ cm

Substituting the values we know into the formula for the arc length,

\begin{aligned} \text{Arc length } &=\frac{\theta}{360} \times 2 \times \pi \times r \\\\ &=\frac{40}{360} \times 2 \times \pi \times 7 \\\\ &=4.886921906 \mathrm{~cm} \end{aligned}

5. Find the length of the arc of the sector. Give your answer to 1 decimal place.

16.8 \ cm

4.2 \ cm

8.4 \ cm

12.6 \ cm

Substituting the values we know into the formula for the arc length, remembering \theta =60^{\circ} because this is an equilateral triangle.

\begin{aligned} \text{Arc length } &=\frac{\theta}{360} \times 2 \times \pi \times r \\\\ &=\frac{60}{360} \times 2 \times \pi \times 7 \\\\ &=8.37758041 \mathrm{~cm} \end{aligned}

6. Find the angle of the sector with an arc length of 10 \pi \ cm and a radius of 12 \ cm.

150^{\circ}

210^{\circ}

60^{\circ}

30^{\circ}

Substituting the values we know into the formula for arc length,

\frac{\theta}{360} \times 2 \times \pi \times 12=10 \pi

Solving for the angle \theta,

\theta=150^{\circ}.

### Arc of a circle GCSE questions

1. A sector has a radius of 15 \ cm and an angle of 120^{\circ}. Find the length of the arc in terms of \pi.

(2 marks)

Use of \frac{\theta}{360} \times \pi d .

(1)

10 \pi

(1)

2. The sector has an angle of 50^{\circ} and an arc length of 25 \ cm. Find the radius of the sector in centimetres to the nearest millimetre.

(3 marks)

25=\frac{50}{360} \times \pi d (or equivalent).

(1)

Diameter = 57.29…

(1)

(1)

3. A metal sector with radius 6 \ cm and angle 250^{\circ} is used to form an open cone, as shown.

Calculate the radius, r, of the cone formed from the sector.

(4 marks)

Use of \frac{250}{360} \times 12 \pi .

(1)

Arc length = 26.1799…

(1)

\pi r^{2}=26.1799…

(1)

r=2.89 \ cm

(1)

## Learning checklist

You have now learned how to:

• Identify and apply circle definitions and properties, including tangent, arc, sector and segment
• Calculate the length of an arc
• Calculate a property of a circle given the length of an arc

## Still stuck?

Prepare your KS4 students for maths GCSEs success with Third Space Learning. Weekly online one to one GCSE maths revision lessons delivered by expert maths tutors.

Find out more about our GCSE maths revision programme.

x

#### GCSE Maths Papers - November 2022 Topics

Practice paper packs based on the November advanced information for Edexcel 2022 Foundation and Higher exams.

Designed to help your GCSE students revise some of the topics that are likely to come up in November exams.


text:
Leon Bilton

# How to not hate LaTeX

##### Last modified on May 17, 2022

In this article I will collect anecdotes and thoughts on $$\LaTeX$$ and similar software. Although this tool is routinely used to create scientific and technical documents, the complexity of the task combined with certain idiosyncrasies of the software mean that it is notoriously uncomfortable to learn and can seem unwieldy. After several years of this experience, I have come to appreciate $$\LaTeX$$ more, not least for its free and open source nature. I find the word like too strong in this context, but I don't hate $$\LaTeX$$. Maybe these notes can also help others to arrive at this point. The article is structured more or less like an introduction to $$\LaTeX$$, and doesn't assume much background.

Contents:

## In the beginning was the processor

Advocates of $$\LaTeX$$ will admit that it does not provide a "conventional" editing experience. Convention is in this case defined by a so-called market share of desktop operating systems, with one fairly obvious caveat being that Linux—the platform of choice for many $$\LaTeX$$ users—is free. A second caveat is that $$\LaTeX$$ can be used on other operating systems. Nevertheless, the idea of having to compile a document before visually affirming the formatting is alien to users of "word processors". It might even seem like a step backwards. That nice "Export to PDF" button (a recent addition in fact) is so much more in vogue. Unfortunately, just like compiling, it is also a one-way trip. Recipients of your PDF document have no reliable way to convert back, not without opening their wallet at least.

The approach here is simply different. A .tex file is a quite portable and efficient format for both the textual content, and the formatting directives of a document. Among other things, this allows for embedding almost arbitrarily complex mathematical expressions. Because .tex files are entirely encoded in plain text, they are tiny and can be transported with ease. Unfortunately a simple "track changes" tool for these files is sorely lacking, which makes collaboration with people who are unfamiliar with $$\LaTeX$$ difficult. I do not think there is a good reason for this continued omission, nor do I think it is necessarily impossible to encode that kind of information inside the .tex file itself—see CriticMarkup for inspiration.

I digress. You still want to see the formatting. At this point, you would normally be directed to either something called "LyX", which attempts to display limited formatting of .tex documents in real time, or an online service like Overleaf which shows a periodically updated preview of the PDF. I recommend neither. While Overleaf is interesting for its collaboration features and integration with academic publishing, I don't think it is a necessary or even good way to introduce $$\LaTeX$$. The point is that you don't need to care about how it looks, except for when you're either done writing or the compiler is shouting at you in gibberish. Getting over this is the first step to learning how to use this software.

Installing a TeX distribution is usually straightforward. I use TeX Live on Linux, but cannot make recommendations for other platforms. If you are new to open source software, please be prepared for a little troubleshooting. Something something… free lunch.

## TeX, LaTeX, LuaLaTeX, ConTeXt, XeTeX, BibTeX…

Stop. Here's a simple rule. Are you writing a paper? If so, then you have to use $$\LaTeX$$ (actually pdflatex), and there's probably a template somewhere, so in this case do check Overleaf or the publisher's website. Don't spend much time with formatting, they are going to change it anyway. The only other one of those you care about is BibTeX, it's a metadata format for citations. But you probably already know all that.

If not, then I suggest to use LuaLaTex. It handles fonts better and is generally not a dinosaur (you can actually type real Unicode characters like ä in the .tex file). It comes with TeX Live by default, so usually you can just run lualatex instead of latex. The Lua- part means that it supports executing Lua code in the .tex file. Not directly mind you, you need to use the \directlua command. I've never used it yet, the UTF-8 support is what sold it for me. No, it is not directly compatible with old $$\LaTeX$$. But you only ever need the latter for papers, in which case the publisher's template will probably be the main restriction anyway.

ConTeXt is a different attempt at a change, but it seems abandoned. Some of it got incorporated in to LuaTex which powers LuaLaTex. XeTeX is… something? It also does UTF-8, but differently? I don't really know, but it's older than LuaTeX. LuaLaTeX also allows the use of some new or improved packages, like polyglossia for multilingual documents or unicode-math to interpret Unicode input in mathematics. No more \varepsilons taking up space and causing verbose, hard to read code for equations. Now you can input ε directly (using an appropriately configured editor), which also works in $$\KaTeX$$ on the web, by the way.

## KOMA-Script

A lot of examples start with \documentclass{article}. But the built-in document classes are a nuisance. For example, you have to load a separate package just to change margins. Luckily some German folks figured that it was a silly idea, and if people are going to create a package for it anyway, why not add at least basic geometry options to the document class. A good starting point is

\documentclass[twocolumn,12pt,DIV=14]{scrarticle}

Remove the twocolumn option to use a single column format, but bear in mind that two column documents are easier to read, and use less paper. The DIV parameter controls the text area. It works like you'd expect, but for recommended font size and DIV combinations, scroll to Table 2.1 in the KOMA-Script manual (~3M PDF download).

Want more space between the two columns?

\setlength{\columnsep}{1cm}

The KOMA-Script classes are also nice once you become more familiar with $$\LaTeX$$, since they expose a lot more options than the default classes. So far, I haven't had any trouble with using them from LuaLaTex. They save a lot of hassle, and eliminate the need for many band aid packages. If you need to change something, first check in the KOMA-Script manual. The option might already be there.

## The preamble

The preamble of a $$\LaTeX$$ document is what goes before the actual document environment, i.e. before \begin{document}, and consists mostly of loading required components and setting global options. A lot of examples out there have long and arcane preambles that would cause even the bravest novice pause. The worst offenders usually include comments that, rather than explaining anything, warn against the hubris of attempting to modify any of the intricate incantations. It should be obvious that this is not a viable strategy, unless you are prepared to house a resident guru.

Instead, it would be better to not put things in the preamble, until you understand what they are doing and can decide that they are really necessary. This is a good way to think about software dependencies in general. I refer the interested reader to the essay by Russ Cox, "Our Software Dependency Problem".

The LuaLaTex manual itself recommends polyglossia for multilingual documents. This means that it provides things like the correct settings for dates and hyphenation patterns (telling LuaLaTex where it should break words to justify text). It also appropriately changes the result of commands like \section, \chapter, etc.

I have already mentioned unicode-math, but an essential package for mathematics is amsmath. This provides things like aligned or line-wrapped equation environments and proper typography of bold math symbols. There is also \text for putting upright text inside an equation.

As a third essential package, I recommend microtype. This provides well researched inter-word spacing adjustments and other fine-tunings, which I generally don't care about but luckily there is no need to set anything manually. Just load the package, and maybe use the final option to ensure that the typographical corrections are applied even in draft mode.

\usepackage[final]{microtype}

This shouldn't be very noticeable, but it can shave off a page or two in long documents.

You'll find that a few more things are needed depending on context, like graphicx for images and hyperref if you want clickable links. The csquotes package can be used for non-English quotation marks (for correct English quotation marks use the unicode characters, i.e. “foo” or ‘foo’). Citations and bibliographies are another story, unfortunately.

Here's an example of a simple but complete preamble (although without bibliography support):

\documentclass[twocolumn,12pt,DIV=15,abstract=true]{scrarticle}
\setlength{\columnsep}{1cm}
\usepackage{polyglossia}
\setdefaultlanguage[variant=british]{english}
\usepackage{libertinus}
\usepackage[final]{microtype}

\usepackage{unicode-math}  % Better math font config, and unicode glyph math input
\usepackage{amsmath}  % Bold symbols, \DeclareMathOperator and $\text{…}$
\usepackage{siunitx}  % Physical units, number ranges, scientific and complex number notation macros

\usepackage{graphicx}  % For \includegraphics
\usepackage{subcaption}  % For the subcaptionblock environment (replaces minipage for subfigures)
\addtokomafont{caption}{\small}  % Use small font for caption text.
\setkomafont{captionlabel}{\bfseries}  % Use bold font for caption label.
\setcapindent{0pt}  % Don't use a hanging indent for multiline captions.

\usepackage{xcolor}  % Without options just provides the ! mechanism, e.g. black!60 for grayscale
\usepackage{xurl}  % Better URL line breaks.
\usepackage[colorlinks=true,allcolors=black!60]{hyperref}  % Clickable references and URLs.
\usepackage[noabbrev]{cleveref}  % References figures/tables/equations properly.

## Intermission: a rant about LaTeX syntax

Ah, the $$\LaTeX$$ syntax. It's just not good. The verbosity is overwhelming, and many commands are not very mnemonic. There are already some hints of this in the provided preamble example. Why is the column separation set with \setlength, but the caption font is changed with \addtokomafont? And then \setkomafont is used for the caption label setting! Unfortunately, this is an example of relatively nice $$\LaTeX$$ syntax. Matrices are written like:

\begin{bmatrix}
a & b \\
c & d
\end{bmatrix}

That's just horrible. Although using unicode to enter symbols alleviates some of the pain, the syntax is a wart that will remain and continue to discourage new users.

It's not just the verbosity or number of commands either. The mathematical "grammar" is all over the place. Some commands are named for their appearance, like \bigotimes or \circ (\deg also exists but just prints the upright letters 'deg'), but others are named for their mathematical function, like \subset or \to. To type the paragraph symbol § we use \S, but for ⋁ we use \bigvee spelled out (there is no \V command that I know of). Combinations of characters like <= are not interpreted, leading to a wonderful plethora of commands like \ll, \leq (or just \le if you're in a rush), \neq, \simeq,… Because the backslash is commonly used for escaping non-printed characters, we also have to do \mid instead of, say, \|.

My personal favourite is \ni. Can you guess what it does?

Finally, and perhaps most importantly, the syntax does not preserve mathematical semantics. In other words, all of the effort spent learning this "language" will be useful for little more than pretty-printing mathematics. It would be cool if equations could be pasted into a computer algebra system after copying them directly out of a $$\LaTeX$$ document, but that won't happen. At least, not without carrying an extremely complex parser around and specifying a bunch of additional disambiguations. For a better attempt at mathematical syntax, see Mathematica. Here is a transcript of a 2010 talk by the man himself. There are also a few interesting open source projects in this area, like Cortex and Metatheory which I may look into at some point.

## Single columns in two-column documents

I have already recommended the use of two-column documents. This will quickly lead to the painful realisation that equations are not line-wrapped automatically. This is a case where $$\LaTeX$$ is simply too dumb. That is one of the reasons for loading the amsmath package, which provides the multiline, gather, split and align environments. They are all useful in different contexts, and starred versions usually suppress the equation number as usual. The folks at Overleaf provide a quick overview.

Luckily, it is possible to get full-width figures and tables, by simply using the figure* and table* environments. However, arbitrary blocks that span the full width of the page require an obscene syntax:

\documentclass{twocolumn,12pt,DIV=14,abstract=true]{scrarticle}
\usepackage[final]{microtype}
\usepackage{amsmath}

\title{Single column abstract in a two-column document}

\begin{document}

\twocolumn[%
\maketitle
\begin{@twocolumnfalse}
\noindent\hspace{0.1\linewidth}\begin{minipage}{0.8\textwidth}

\begin{abstract}
Hey look! We can finally type the abstract in here.
It will appear centered and spread across 80\% of the page.
\end{abstract}

\end{minipage}
\bigskip
\bigskip
\end{@twocolumnfalse}]

\end{document}

So to initialise the single-column box we use \twocolumn. A pretty good start, right? Oh, and here you need to use square brackets, because we are putting stuff into the "optional" argument of that command. I don't care to find out why.

The \maketitle command must be put inside the block, otherwise a page break appears after the title. Then, we use the actual @twocolumnfalse environment. The leading @ is required, it's just how things are.

You might have imagined that for centering the abstract, a command like \centering would be quite useful. This command does in fact exist, but when it is used inside \twocolumn it doesn't work as you'd expect. Instead, we shift the box manually by 0.1\linewidth (10% of the page width), taking care to enforce that the box itself is not indented like a paragraph, and then set up a minipage spanning 80% (note the use of \textwidth) of the page.

The \bigskips at the end make sure there is appropriate spacing between the box and the following content, which is somehow not the default.

Overall, it only took a triple-nested environment and five different commands. Success?

There is one more caveat to this approach. Because of how $$\LaTeX$$ handles “floating” boxes like this, any so-called “fragile” commands, like \cite for example, will need to be enclosed in a \protect command to avoid mysterious and entirely unhelpful errors like

Argument of \blx@citeargs@i has an extra }.

when trying to do \parencite[e.g.][]{key} inside the box. Instead use \protect{\parencite[e.g.][]{key}}. The same applies for citations inside figure captions.

## Bibliographies, citations and references

References to figures or tables are generally fairly easy, and the cleveref package simply remove the need to type 'table' or 'figure', instead we can just do

This is a reference to~\cref{<label>}.

where <label> is the label given to the figure or table via the \label command. The tilde ~ in front is used to tell the software to never break the line at that spot. To capitalise the word used by cleveref, e.g. at the start of a sentence, use \Cref instead.

Citations and bibliographic entries are more involved to set up, but they are nicely automated. The most annoying part is acquiring the reference metadata, which must be in BibTeX format. The BibTeX format has undergone various revisions and extensions, which can also complicate things. These days, journals will often provide a BibTeX snippet for most papers, if you can find it. It is also possible to search for the paper directly on CrossRef, and they usually provide higher quality BibTeX (the snippets from journals can be incomplete or garbled). That is to be expected, because these are the people who register Digital Object Identifiers (DOIs) for most academic papers out there and some other content as well. I have written a small python script to fetch BibTeX snippets automatically, given the DOI of a paper. It doesn't always work, and not all papers even have a DOI. Be prepared to manually fill in or correct missing BibTeX occasionally.

The BibTeX snippets are collected in a file ending with .bib. There are various ways to load this file into $$\LaTeX$$. Whenever possible, I use the newer biber engine, and add something like this to the preamble:

\usepackage[backend=biber,bibencoding=utf8,style=apa]{biblatex}
\renewcommand*{\bibfont}{\normalfont\small}
\addbibresource{\jobname.bib}

The command \jobname resolves to the base name of the .tex file, so I just name my .bib file identically. Citations are included in the text by using the \cite or \parencite commands.

Now, because $$\LaTeX$$ is pretty stupid, it doesn't handle citations properly on its own. Not only do we need to load the biblatex package (or similar, e.g. the legacy natbib), but we also need to run a different command on the output of the first lualatex run. For me, this is biber, as indicated by the backend option in the snippet.

In fact, it doesn't stop there. The biber tool only inserts citation text in some kind of naïve way or whatever, so we need to run lualatex again afterwards. No, no we're not done. To be completely sure that the formatting is readjusted after inserting citations, bibliography entries and possibly back-references, it is actually necessary to run lualatex twice afterwards.

Wow, that's pretty stupid. It's so stupid that a variety of build tools have sprung up to properly automate the document generation. That's right build tools, you are really just a programmer now. There are really two options: find an editor that has a $$\LaTeX$$ plugin and can properly automate this absurd workflow, or embrace your inevitable fate…

src := $(wildcard *.tex) bib :=$(wildcard *.bib)
all: $(src:%.tex=%.pdf) %.pdf:$(src) out/%.bbl
@ 1>/dev/null lualatex --output-directory=out --interaction=batchmode --halt-on-error --file-line-error $< \ || echo "LuaLaTeX run 2/3 failed, check log files in out/ for details" @ lualatex --output-directory=out --halt-on-error --file-line-error$<|grep -v '/usr/share/'
@ mv out/$@$@

out/%.bbl: out/%.bcf
@ biber --input-directory=out --output-directory=out $(<F) out/%.bcf:$(src) $(bib) | out @ 1>/dev/null lualatex --output-directory=out --interaction=batchmode --halt-on-error --file-line-error$< \
|| echo "LuaLaTeX run 1/3 failed, check log files in out/ for details"

out:
@ mkdir -p out

This is a Makefile example for LuaLaTeX + Biber. You thought that $$\LaTeX$$ syntax was bad, eh? I use it by storing the .tex and .bib files in the same folder, and just running make to compile the document. It's a Linux thing, I won't go over it here. Check GNU Make if you're interested. If you use a Linux platform, it should normally be installed already.

## Overfull and underfull, isn't it just wonderfull

With any $$\LaTeX$$ setup, you will eventually encounter the infamous warnings about 'overfull' or 'underfull' boxes. They are quite terse and mysterious:

Underfull \hbox (badness 1320) in paragraph at lines 691--691

Hmm, my .tex file doesn't even have 691 lines. That's right, don't waste time looking at the line number, it refers to some intermediate representation of the content. Instead, check the lines immediatelly following this in the command output (which is far to messy, another gripe). Usually these warnings are caused by

• using long, techincal terms that the software doesn't know how to hyphenate, which prevents proper line breaks, or

• using inline equations that put the inter-word spacing out of whack, or

• not wrapping an equation with the amsmath environments, or

• doing something that the software doesn't like such as setting a block across the whole page in a two-column document

Avoiding these warnings is not critical and sometimes not practical. If there are some clues after the error, glance over the PDF and check for weird spacing issues. To add hyphenation hints to a long word, use \-, e.g. whats\-amajig. If word spacing is too large because of an inline equation, it is sometimes possible to fix the warning by adding e.g. \hspace{4pt} before and after the equation. For the two-column issue, I don't know of any way to prevent the warnings.

## Fin

That's all for now. Go away and compile your document. Yes I know, the compiler is too slow. Parallel $$\LaTeX$$ compilation, you ask? Hah, check back in a few decades, if you're extremely lucky we'll have multi-threading…