text:
# Age prediction using machine learning

### Overview

question Questions
• How to use machine learning to create predictive models from biological datasets (RNA-seq and DNA methylation)?

objectives Objectives
• Learn ageing biomarkers from RNA-seq and DNA methylation datasets

• Apply regression based machine learning algorithms

• Learn feature selection and hyperparameter optimisation

requirements Requirements

time Time estimation: 2 hours

Supporting Materials

# Introduction

Machine Learning is used to create predictive models by learning features from datasets. In the studies Jason G. Fleischer et al. 2018 and Jana Naue et al. 2017, biomarkers are examined to predict the chronological age of humans by analysing the RNA-seq gene expression levels and DNA methylation pattern respectively. Different machine learning algorithms are used in these studies to select specific biomarkers to make age prediction. The RNA-seq gene expression (FPKM) dataset is generated using fibroblast cell lines of humans. The skin fibroblasts cells keep damages that happen with age. Epigenomic and phenotypic changes which are age-dependent are also contained in these cells. Within each individual, DNA methylation changes with age. This knowledge is used to select useful biomarkers from DNA methylation dataset. The CpGs sites with the highest correlation to age are selected as the biomarkers/features. In both these studies, specific biomarkers are analysed by machine learning algorithms to create an age prediction model.

The datasets from these studies contain features (present as columns). The last column in both the datasets refers to age which is used as labels/targets. Since the targets are real numbers, the machine learning task becomes regression. Using these features and targets, a model is learned using machine learning (regression) which learns a mapping between these features and targets. Using machine learning tools in Galaxy, we can achieve comparable prediction scores as achieved by these analyses. In this tutorial, we will apply a couple of (scikit-learn) machine learning tools to RNA-seq and DNA methylation datasets to predict the chronological age of humans. This tutorial is divided into two parts - one with RNA-seq dataset and another with DNA methylation dataset.

### comment Comment

In cross-validation, a dataset is divided into k equal parts (k should be at least 2). One part is used as the validation set to evaluate the performance of any machine learning algorithm and the remaining parts are used for learning/training.

#### R2 (coefficient of determination)

In both the parts, learning on datasets is done using cross-validation and R2 scoring metric is used to evaluate the performance of the trained model. The closer it is to 1.0, the better it is. If it is negative, then the trained model is not good. To infer how its values exhibit model performance, we can compare the figures 1 and 2. In both the plots, the true and predicted targets are plotted in a scatter plot. For a good model, most of the points should lie along the x = y line as the true and predicted targets are close to each other. In figure 1, we can see that the points are scattered and do not show any pattern. Therefore, the R2 score is -0.06. But, figure 2 shows a better pattern as most of the points lie along the line and the R2 score is almost 1.0. For RNA-seq dataset, we will compute the cross-validated R2 score using the training set and for DNA methylation dataset, we will compute the R2 score for the test set.

### Agenda

In this tutorial, we will cover:

1. Analyze RNA-seq dataset
2. Analyze DNA methylation dataset

# Analyze RNA-seq dataset

The RNA-seq dataset is collected from fibroblast cell lines belonging to 133 healthy patients with age ranging from 1 to 94 years. Skin fibroblasts cells dataset is conducive for age prediction related studies because of multiple reasons - these skin fibroblasts hold age-related damages, contain phenotypic, epigenomic and transcriptomic changes which are age-dependent and it is easy to collect this dataset using non-invasive techniques. Using this dataset, we perform an exhaustive search (also known as grid search) for finding the best features and then apply ElasticNet regressor with 5-fold cross-validation. The R2 score achieved is comparable to the predictions found in the study Jason G. Fleischer et al. 2018.

### details 5-fold cross-validation

It is a model validation technique which estimates the performance of a predictive model on an unseen data. A dataset is divided into 5 folds and these folds are categorised into training and validation sets. The idea of cross-validation is shown in figure 3. The complete dataset is divided into 5 equal parts. 80% of the dataset is used for training and the remaining 20% is used for validating the performance of training. This is done for 5 folds/iterations, each time the validation set (20% of the dataset) is different. In all five folds, the complete dataset is used for training and validation. The final validation performance is averaged over 5 folds.

## Get dataset

We proceed to the analysis by uploading the RNA-seq dataset. The dataset has 133 rows corresponding to humans and over 27,000 columns specifying genes. The last column age contains the chronological age of humans. For the validation set, this age column is predicted and R2 metric is computed.

### hands_on Hands-on: Data upload

1. Create a new history for this tutorial

### tip Tip: Creating a new history

Click the new-history icon at the top of the history panel

If the new-history is missing:

1. Click on the galaxy-gear icon (History options) on the top of the history panel
2. Select the option Create New from the menu
2. Import the files from Zenodo

https://zenodo.org/record/2545213/files/training_data_normal.tsv

• Copy the link location
• Open the Galaxy Upload Manager (galaxy-upload on the top-right of the tool panel)

• Select Paste/Fetch Data
• Paste the link into the text field

• Press Start

• Close the window

By default, Galaxy uses the URL as the name, so rename the files with a more useful name.

3. Rename the dataset to training_data_normal.

### tip Tip: Renaming a dataset

• Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
• In the central panel, change the Name field
• Click the Save button
4. Check that the datatype is tabular

### tip Tip: Changing the datatype

• Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
• In the central panel, click on the galaxy-chart-select-data Datatypes tab on the top
• Select datatypes
• Click the Change datatype button

### comment Comment

In the training_data_normal dataset, you might have noticed that it showed 134 rows instead of 133 rows. The first row contains the header information i.e. the description of each column.

## Create data processing pipeline

We can see that this RNA-seq dataset is high-dimensional. There are over 27,000 columns/features. Generally, not all the features in the dataset are useful for prediction. We need only those features which increase the predictive ability of the model. To filter these features, we perform feature selection and retain only those which are useful. To do that, we use SelectKBest module. This approach involves extracting those features which are most correlated to the target (age in our dataset). F-regression is used for the extraction of features. Moreover, we are not sure of how many of these features we will need. To find the right number of features, we do a hyperparameter search (finds the best combination of values of different parameters). It works by setting a different number of features and find out the number for which the accuracy is the best among all the numbers. To wrap this feature selector with a regressor, we will use the Pipeline builder tool. This tool creates a sequential flow of algorithms to execute on datasets. It does not take any dataset as input. Rather, it is used as an input to the Hyperparameter search tool (explained in the following step). We will use ElasticNet as a regressor which creates an age prediction model. It is a linear regressor with l1 (also called lasso) and l2 (also called ridge) as regularisers. Regularisation is a technique used in machine learning to prevent overfitting. Overfitting happens when a machine learning algorithm starts memorising dataset it is trained upon instead of learning general features. The consequence of overfitting is that the accuracy on the training set is good but on the unseen set (test set) is not good which happens because the algorithm has not learned general features from the dataset. To prevent overfitting, regularisers like l1 and l2 are used. L1 is a linear term added to the error function of a machine learning algorithm and l2 adds a squared term to the error function. More details about l1 and l2 can found here.

### hands_on Hands-on: Create pipeline

1. Pipeline builder tool with the following parameters:
• In “1: Pre-processing step”:
• “Choose the type of transformation”: Feature Selection
• “Select a feature selection algorithm”: SelectKBest - Select features according to the k highest scores
• “Select a score function”: f_regression - Univariate linear regression tests
• In “Final Estimator”:
• “Choose the module that contains target estimator”: sklearn.linear_model
• “Choose estimator class”: ElasticNet
• In “Output the final estimator instead?”: Pipeline

## Optimise hyperparameters

In any machine learning algorithm, there are many parameters (hyperparameters). We are not sure which values of these parameters will give an optimal prediction. The default values given for these parameters may not be optimal for different datasets. To find the best combination of the values of different parameters for a dataset, hyperparameter optimisation is performed. There are different techniques to optimise the hyperparameters of any algorithm given a dataset:

For our analyses, we will use the grid search approach. It is an exhaustive search which tries out all the combinations of different hyperparameters and ranks these combinations based on a scoring metric. In the random search, the values of a parameter are selected randomly from a given range and the best one is found. Grid search works well for parameters taking categorical as well as numerical values but for the random search, it becomes difficult for parameters which take categorical values.

In the Pipeline builder tool, we added two steps - preprocessing (feature selection) and an estimator (regressor). There are different hyperparameters for these two steps and their best combination should be found out. We will perform grid search to estimate the best values for these parameters: k (number of features), normalize (subtract the mean and divide by the l2-norm of the dataset) and alpha (a constant which is multiplied to the regularisation term). For each parameter, we need to specify a set of values to choose from:

• k: [5880, 5890, 5895, 5900]

These values of k are chosen to get the best accuracy. We can choose any number (integers) between 1 and 27,000 (maximum number of features in the dataset). We will use only these values (shown above) for k as the accuracy remains the best around these numbers. But, it may vary for a different RNA-seq dataset. That’s the reason why we perform hyperparameter search to find the best values of parameters for any dataset.

• normalize: [True, False]

The default value of normalize is False. We will check both, True and False.

• alpha: [0.00001, 0.0001, 0.001]

The parameter alpha takes a positive real number and its default value is 1.0.

For these three parameters, we have 24 different combinations (4 x 2 x 3) of values and we will verify the performance of each combination. The parameter k is used for feature selection and parameters normalize and alpha are used for regressor. There are many more hyperparameters of ElasticNet regressor which are explained in the official documentation of scikit-learn. But, the combination of the above three parameters already gives a comparable accuracy published in the study Jason G. Fleischer et al. 2018. Therefore, we will stick to these parameters.

### comment Comment

It is advisable to tune all the parameters of a machine learning algorithm for a dataset if no prior information is available about the subset of parameters which works best for the dataset.

### comment Comment

These parameters have the same description and values in the second part of the tutorial where we will again use the Hyperparameter search tool.

### Extract hyperparameters

Before searching for the best values of hyperparameters, we require a tool to extract the list of hyperparameters of data preprocessors and estimators. To achieve it, we will use the Estimator attributes tool. This tool creates a tabular file with a list of all the different hyperparameters of preprocessors and estimators. This tabular file will be used in the Hyperparameter search tool to populate the list of hyperparameters with their respective values.

### hands_on Hands-on: Estimator attributes

1. Estimator attributes tool with the following parameters:
• param-files “Choose the dataset containing estimator/pipeline object”: pipeline builder file (output of Pipeline builder tool)
• “Select an attribute retrieval type”: Estimator - get_params()

### Search for the best values of hyperparameters

After extracting the parameter names from the Pipeline builder file using Estimator attributes tool, we will use the Hyperparameter search tool to find the best values for each hyperparameter. These values will lead us to create the best model based on the search space chosen for each hyperparameter.

1. Hyperparameter search tool with the following parameters:
• “Select a model selection search scheme”: GridSearchCV - Exhaustive search over specified parameter values for an estimator
• param-files “Choose the dataset containing pipeline/estimator object”: zipped file (output of Pipeline builder tool)
• In “Search parameters Builder”:
• param-files “Choose the dataset containing parameter names”: tabular file (output of Estimator attributes tool)
• In “Parameter settings for search”:
• param-repeat “1: Parameter settings for search”
• “Choose a parameter name (with current value)”: selectkbest__k: 10
• “Search list”: [5880, 5890, 5895, 5900]
• param-repeat “2: Parameter settings for search”
• “Choose a parameter name (with current value)”: elasticnet__normalize: False
• “Search list”: [True, False]
• param-repeat “3: Parameter settings for search”
• “Choose a parameter name (with current value)”: elasticnet__alpha: 1.0
• “Search list”: [0.00001, 0.0001, 0.001]
• In “Advanced Options for SearchCV”:
• “Select the primary metric (scoring)”: Regression -- 'r2'

A scoring metric can be set. In this tutorial, we use Regression -- 'r2'

• “Select the cv splitter”: KFold

There are different ways to split the dataset into training and validation sets. In our tutorial, we will use KFold which splits the dataset into K consecutive parts. It is used for cross-validation. It is set to 5 using another parameter n_splits.

• “n_splits”: 5
• “Whether to shuffle data before splitting”: Yes
• “Random seed number”: 3111696

It is set to an integer and used to retain the randomness/accuracy when “Whether to shuffle data before splitting” is Yes across successive experiments.

• “Raise fit error”: No

While setting different values for a parameter during hyperparameter search, it can happen that wrong values are set which may generate exceptions. To avoid stopping the execution of a regressor, it is set to No which means even if a wrong parameter value is encountered, the regressor does not stop running and skips that value.

• “Select input type”: tabular data
• param-files “Training samples dataset”: training_data_normal tabular file
• “Does the dataset contain header”: Yes
• “Choose how to select data by column”: All columns BUT by column header name(s)
• “Type header name(s)”: age
• param-files “Dataset containing class labels or target values”: training_data_normal tabular file
• “Does the dataset contain header”: Yes
• “Choose how to select data by column”: Select columns by column header name(s)
• “Type header name(s)”: age

The tool returns two outputs, one of which is a table with numerical results. Please inspect it carefully: the rank_test_score column shows the ranking of different combinations based on the values in the mean_test_score column.

### question Questions

1. What is the best mean_test_score value estimated by the Hyperparameter search tool?
2. Which combination of parameters gives the best result?
3. How many possible combinations of parameters the Hyperparameter searchh tool estimated?

### solution Solution

1. 0.73 (it is close to the best R2 score (0.81) achieved by a customised ensemble algorithm explained in Jason G. Fleischer et al. 2018)
2. alpha: 0.001, normalize: True, k: 5880
3. 24 (it is equal to the number of rows in the tabular output of Hyperparameter search tool)

## Create parallel coordinates plot

We will visualize the tabular output of Hyperparameter search tool from the previous step using the parallel coordinates plot of tabular data. There are multiple columns in the tabular output, but we will focus on only a few of them.

### hands_on Hands-on: Create parallel coordinates plot

1. Parallel coordinates plot tool with the following parameters:
• param-files “Select data file”: tabular file (output of Hyperparameter search tool)
• “Select the columns for dimensions”: c4, c5, c6
• “Select a column containing the values for coloring”: c3

The output plot has the following legend: the colour-coding is based on the mean_test_score (c3) column. You can follow the line leading to the score along every column with parameters’ settings. The columns c4, c5 and c6 are the parameters we chose and c3 is the accuracy column present in the tabular output of Hyperparameter search tool.

### question Questions

What are the worst performing combinations of hyperparameters (name four)?

### solution Solution

Worst four:

• alpha: 0.00001, normalize: False, k: 5880
• alpha: 0.00001, normalize: False, k: 5890
• alpha: 0.00001, normalize: False, k: 5895
• alpha: 0.00001, normalize: False, k: 5900

## Summary

Figure 4 shows that we achieved an R2 score of 0.73 (last column) with 5-fold cross-validation on the training set. In the study Jason G. Fleischer et al. 2018 as well, a similar R2 score is mentioned for linear regressors (Linear regression and ElasticNet). Moreover, the study also included a customised ensemble regressor which achieved better performance (R2 = 0.81). However, our analysis showcases the use of machine learning tools in Galaxy to reproduce the results published in the paper.

# Analyze DNA methylation dataset

In the second part of the analysis, we will use the DNA methylation dataset to predict chronological age. One important reason to choose this dataset for age prediction task is that DNA methylation changes with age and this change occurs at specific CpG sites in humans. Whole blood samples are collected from humans with their ages falling in the range 18-69 and the best age-correlated CpG sites in the genome are chosen as features. The dataset is divided into two parts - training and test sets. The training set is used to train a regressor and the test set is used to evaluate the performance of the trained model using the R2 scoring metric. 5-fold cross-validation is used for training.

## Get train and test datasets

We proceed with the analysis by uploading new datasets. You might want to create a new history first. The training set contains 208 rows corresponding to humans and 13 features (age-correlated CpG sites in DNA methylation dataset). The last column is age. The test set contains 104 rows and the same number of features as the training set. The age column in the test set is predicted after training on the training set. Another dataset test_rows_labels contains the true age values of the test set which is used to compute R2 scores between true and predicted age.

### hands_on Hands-on: Data upload

1. Create a new history for this tutorial
2. Import the files from Zenodo

https://zenodo.org/record/2545213/files/train_rows.csv
https://zenodo.org/record/2545213/files/test_rows_labels.csv
https://zenodo.org/record/2545213/files/test_rows.csv

• Copy the link location
• Open the Galaxy Upload Manager (galaxy-upload on the top-right of the tool panel)

• Select Paste/Fetch Data
• Paste the link into the text field

• Press Start

• Close the window

By default, Galaxy uses the URL as the name, so rename the files with a more useful name.

3. Rename the datasets as train_rows, test_rows_labels and test_rows respectively.

### tip Tip: Renaming a dataset

• Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
• In the central panel, change the Name field
• Click the Save button
4. Check that the datatype of all the three datasets is tabular.

### tip Tip: Changing the datatype

• Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
• In the central panel, click on the galaxy-chart-select-data Datatypes tab on the top
• Select datatypes
• Click the Change datatype button

The train_rows contains a column Age which is the label or target. We will evaluate our model on test_rows and compare the predicted age with the true age in test_rows_labels

## Create data processing pipeline

We will create a pipeline with Pipeline builder tool but this time, we just specify the regressor. Jana Naue et al. 2017 has used Random Forest as the regressor and we can conclude from this study that the ensemble-based regressor works well on this DNA methylation dataset. Therefore, we will use Gradient boosting which is an ensemble-based regressor because it uses multiple decision tree regressors internally and predicts by taking the collective performances of the predictions (by multiple decision trees). It has a good predictive power and is robust to the outliers. It creates an ensemble of weak learners (decision trees) and iteratively minimises error. One disadvantage which comes from its basic principle of boosting is that it cannot be parallelised. The Pipeline builder tool will wrap this regressor and return a zipped file. We will use this zipped file with Estimator attributes tool set the search space of hyperparameters.

### hands_on Hands-on: Create pipeline

1. Pipeline builder tool with the following parameters:
• In “Final Estimator”:
• “Choose the module that contains target estimator”: sklearn.ensemble
• “Choose estimator class”: GradientBoostingRegressor
• In “Output the final estimator instead?”: Final Estimator

We choose Final Estimator as we have only the estimator and no preprocessor and need the parameters of only the estimator.

### comment Comment

Ensemble method uses multiple learning models internally for better predictions.

## Optimise hyperparameters

For this analysis as well, we will use the Hyperparameter search tool to estimate the best values of parameters for the given dataset. We use only one parameter n_estimators of Gradient boosting regressor for this task. This parameter specifies the number of boosting stages the learning process has to go through. The default value of n_estimators for this regressor is 100. But, we are not sure if this gives the best accuracy. Therefore, it is important to set this parameter to different values to find the optimal one. We choose some values which are less than 100 and a few more than 100. The hyperparameter search will look for the optimal number of estimators and gives the best-trained model as one of the outputs. This model is used in the next step to predict age in the test dataset.

### Extract hyperparameters

We will use the Estimator attributes tool to get a list of different hyperparameters of the estimator (including n_estimators). This tool creates a tabular file with a list of all the different hyperparameters of the preprocessors and estimators. This tabular file will be used in the Hyperparameter search tool to populate the list of hyperparameters with their respective (default) values.

### hands_on Hands-on: Estimator attributes

1. Estimator attributes tool with the following parameters:
• param-files “Choose the dataset containing estimator/pipeline object”: final estimator builder file (output of Pipeline builder tool)
• “Select an attribute retrieval type”: Estimator - get_params()

### Search for the best values of hyperparameters

After extracting the parameter names from the Pipeline builder file, we will use the Hyperparameter search tool to find the best values for each hyperparameter. These values will lead us to create the best model based on the search space chosen for each hyperparameter.

### hands_on Hands-on: Hyperparameter search

1. Hyperparameter search tool with the following parameters:
• “Select a model selection search scheme”: GridSearchCV - Exhaustive search over specified parameter values for an estimator
• param-files “Choose the dataset containing pipeline/estimator object”: zipped file (output of Pipeline builder tool)
• In “Search parameters Builder”:
• param-files “Choose the dataset containing parameter names”: tabular file (output of Estimator attributes tool)
• In “Parameter settings for search”:
• param-repeat “1: Parameter settings for search”
• “Choose a parameter name (with current value)”: n_estimators: 100
• “Search list”: [25, 50, 75, 100, 200]
• In “Advanced Options for SearchCV”:
• “Select the primary metric (scoring)”: Regression -- 'r2'

A scoring metric can be set. In this tutorial, we use Regression -- 'r2'

• “Select the cv splitter”: KFold

There are different ways to split the dataset into training and validation sets. In our tutorial, we will use KFold which splits the dataset into K consecutive parts. It is used for cross-validation. It is set to 5 using another parameter n_splits.

• “n_splits”: 5
• “Whether to shuffle data before splitting”: Yes
• “Random seed number”: 3111696

It is set to an integer and used to retain the randomness/accuracy when “Whether to shuffle data before splitting” is Yes across successive experiments.

• “Raise fit error”: No

While setting different values for a parameter during hyperparameter search, it can happen that wrong values are set which may generate exceptions. To avoid stopping the execution of a regressor, it is set to No which means even if a wrong parameter value is encountered, the regressor does not stop running and skips that value.

• “Select input type”: tabular data
• param-files “Training samples dataset”: train_rows tabular file
• “Does the dataset contain header”: Yes
• “Choose how to select data by column”: All columns BUT by column header name(s)
• “Type header name(s)”: Age
• param-files “Dataset containing class labels or target values”: train_rows tabular file
• “Does the dataset contain header”: Yes
• “Choose how to select data by column”: Select columns by column header name(s)
• “Type header name(s)”: Age

### question Question

What is the optimal number of estimators for the given dataset?

Hint: Please look at the mean_test_score column in the tabular result from the Hyperparameter search tool.

### solution Solution

75 (Even though the default value of the number of estimators for Gradient boosting regressor is 100, 75 gives the best accuracy. That’s why it is important to perform hyperparameter search to tune these parameters for any dataset). 50 estimators also give almost the same accuracy.

## Predict age

Using the Hyperparameter search tool, we found the best model based on the training data. Now, we will predict age in the test dataset using this model in order to see if the model has learned important features which can generalise on a new dataset. The test dataset (test_rows) contains the same number of features but does not contain the age column. This is predicted using the trained model.

### hands_on Hands-on: Predict age

1. Ensemble methods for classification and regression tool with the following parameters:
• “Select a Classification Task”: Load a model and predict
• param-files “Models”: zipped file (output of Hyperparameter search tool)
• param-files “Data (tabular)”: test_rows tabular file
• “Does the dataset contain header”: Yes

## Create plots

In the previous step, we generated predictions for the test dataset. We have one more dataset (test_rows_labels) which contains the true age values of the test set. Using the true and predicted values of age in the test set, we will verify the performance by analysing the plots.

### hands_on Hands-on: Create regression plots

1. Plot actual vs predicted curves and residual plots of tabular data tool with the following parameters:
• param-files “Select input data file”: test_rows_labels tabular file
• param-files “Select predicted data file”: tabular file (output of Ensemble methods for classification and regression tool)

The tool outputs three HTML files with the interactive plots.

### question Question

Inspect the plots. What can you say about the predictions?

### solution Solution

Figures 6 and 8 show that the prediction is good because the predicted age lies close to the true age.

We can see in the scatter plot (figure 6) that most of the points lie along the x=y curve. It means that the true and predicted ages are close to each other. The root mean square error (RMSE) is 3.76 and the R2 score is 0.94.

The residual plot shown in figure 7 is generated to see if there is any visible pattern between residual (predicted age - true age) and predicted age. For a good model, there should not be any visible pattern with the plotted points.

## Summary

Figure 6 shows that we achieved an R2 score of 0.94 and root mean square score of 3.76 for the test set using Gradient boosting regressor. In the study Jana Naue et al. 2017 as well, a similar root mean square score (3.93) is mentioned using random forest regressor. The root mean square score shows the difference in the true and predicted age of humans. The R2 score (0.94) is close to the best achievable score of 1.0 which shows that the trained model is good. Overall, the second part of the analysis also shows that using the machine learning tools in Galaxy, we can achieve state-of-the-art predictions mentioned in the recent scientific studies.

# Conclusion

In our tutorial, we could reproduce the results using machine learning tools in Galaxy as achieved in these scientific studies - Jason G. Fleischer et al. 2018 and Jana Naue et al. 2017. We also learned how to work with high-dimensional datasets, perform hyperparameter search and cross-validation. Further, we can reuse the trained models to make predictions on a new dataset provided that this new dataset has the same features. There are numerous other machine learning algorithms available in Galaxy which can also be tried out on these datasets to verify whether the accuracy can be improved.

### keypoints Key points

• Various machine learning algorithms should be used to find the best ones

• For each machine learning algorithm, it hyperparameters should be optimised based on the dataset

• Feature selection should be done for high-dimensional datasets


text:
# Currying and Macros

Consider the following example Racket code:1

(define (gt arr ptr) (list arr (add1 ptr)))
(define (eq arr ptr) (list arr ptr))
(define (lt arr ptr) (list arr (sub1 ptr)))

This is getting a bit repetitive, so our instinct as lisp programmers might be to write a macro that captures the repetition. The only things that change across definitions are the function name, and the thing that it does to ptr, which we can call NAME and OP respectively.

(define-macro (ptr-motion NAME OP)
#'(define (NAME arr ptr) (list arr (OP ptr))))

With this macro, now we can define our above operations in terms of their common functionality, which makes this much simpler code. Also, if we ever want to extend this code (as I did to include eq) it’s better to trust the macro than attempt to rewrite possibly-complicated functionality after not looking at it for a while.

(ptr-motion gt add1)
(ptr-motion eq id)
(ptr-motion lt sub1)

This is actually a case where code-writing macros are a bit overkill, though. All we really need are higher-order functions (so we can pass around OP), and some way to easily create functions from other functions, but with some arguments applied. Because all functions are curried in Haskell, it becomes extremely easy to write this macro as a plain function, and then create the operations we want by partial application to that function.

ptrMotion :: (Int -> Int) -> Vector -> Int -> (Vector, Int)
ptrMotion op arr ptr = (arr, op ptr)

We also lose the NAME argument, since we can just use those names directly.

gt = ptrMotion (+1)
eq = ptrMotion id
lt = ptrMotion (-1)

1. You might need something like this for a functional bf implementation↩

format: plain


text:
My Math Forum (http://mymathforum.com/math-forums.php)
-   Calculus (http://mymathforum.com/calculus/)
-   -   Challenging integral marathon! (http://mymathforum.com/calculus/345077-challenging-integral-marathon.html)

 Joppy October 4th, 2018 06:42 PM

Challenging integral marathon!

The idea is poached from the AOPS forum.

Rules: Solve the integral in the previous post (show your working) and provide the next challenging integral for the next person. This is open to suggestion/refinement, I'm just getting us started..

Problem 1 (by greg1313)

Quote:
 Originally Posted by greg1313 (Post 600352) $$\int\frac{x-1}{x+x^2\log(x)}\,dx$$

Solution 1 (by Skipjack)

Quote:
 Originally Posted by skipjack (Post 600367) For $x$ > 0, $\displaystyle \!\int\! \frac{x - 1}{x + x^2\ln(x)}\,dx = \!\int\!\left(\frac{1 + \ln(x)}{1 + x\ln(x)} - \frac{1}{x}\right)dx = \ln(1 + x\ln(x)) - \ln(x) + \mbox{C}$, where $\mbox{C}$ is a constant.

Problem 2

$$\int_0^\infty x e^{1-x} - \left\lfloor{x}\right\rfloor e^{1 - \left\lfloor{x}\right\rfloor} \, dx$$

 v8archie October 4th, 2018 07:39 PM

\begin{align}\int_0^\infty x e^{1-x} - \left\lfloor{x}\right\rfloor e^{1 - \left\lfloor{x}\right\rfloor} \, \mathrm dx &= \int_0^\infty x e^{1-x} \,\mathrm dx - \left( 1+2e^{-1} + 3e^{-2} + \ldots \right) \\ &= -xe^{1-x}\bigg|_0^\infty + \int_0^\infty e^{1-x} \,\mathrm dx - \left( 1\left(e^{-1}\right)^0+2\left(e^{-1}\right)^1 + 3\left(e^{-1}\right)^2 + \ldots \right) \\ &= 0 - e^{1-x}\bigg|_0^\infty - \left(\frac{\mathrm d}{\mathrm dx} (1-x)^{-1}\right) _{x=e^{-1}} \\ &= e-\left(-(1-x)^{-2}\right) _{x=e^{-1}} \\ &= e + \left( \frac1{1-e^{-1}}\right)^2 \\ &= e +\frac{e^2}{e^2-2e+1} \\ &= \frac{{e^3-2e^2+e}+e^2}{e^2-2e+1} \\ &= \frac{e^3-e^2+e}{e^2-2e+1}
\end{align}

 SDK October 4th, 2018 08:38 PM

Rewrite it as a sum
$\int_0^\infty x^{1-x} - \lfloor x \rfloor e^{1-\lfloor x \rfloor} = \sum_{n=0}^\infty \int_n^{n+1} x^{1-x} \ dx - ne^{1-n}$
Now compute the integral via $\frac{d}{dx}(x+1)e^{1-x} = -xe^{1-x}$ and evaluate on $[n,n+1]$. This simplifies the sum to
$\sum_{n=0}^\infty (2-n)e^{1-n}$
but this is just the series expansion for the derivative of $\frac{e^2}{1-x}$ evaluated at $x = e^{-1}$. Since this series converges uniformly for $|x| < 1$, it follows that $\sum_{n=0}^\infty (2-n)e^{1-n} = -\frac{e^2}{(1-e^{-1})^2}$

 Joppy October 4th, 2018 08:54 PM

Alternate solutions are good but please provide a new problem afterwards so we can keep it going :).

 v8archie October 4th, 2018 09:48 PM

$$\int \frac{\sqrt{x+1} - \sqrt{x-1}}{\sqrt{x+1}+\sqrt{x-1}} \,\mathrm dx$$

Answers without trigonometric or hyperbolic functions.

 SDK October 6th, 2018 07:49 AM

This one is quite good. I still haven't found a solution without trig/complex exponentials. Looking forward to reading the solution when somebody smarter than me figures it out.

 v8archie October 6th, 2018 08:53 AM

I don't mean that you can't use them to get to your answer, but the final answer shouldn't contain them.

 skipjack October 6th, 2018 09:14 AM

Quote:
 Originally Posted by v8archie (Post 600470) Answers without trigonometric or hyperbolic functions.
Any final answer that happens to involve trigonometric or hyperbolic functions explicitly can be rewritten using exponential or logarithmic functions instead, so that's not a significant restriction.

 SDK October 6th, 2018 12:53 PM

Quote:
 Originally Posted by v8archie (Post 600566) I don't mean that you can't use them to get to your answer, but the final answer shouldn't contain them.
I finally have a solution without using any trig/complex exponentials but I don't have a good integral to suggest to replace it so I'll defer answering until I think of something good.

 SDK October 6th, 2018 03:41 PM

Quote:
 Originally Posted by v8archie (Post 600470) $$\int \frac{\sqrt{x+1} - \sqrt{x-1}}{\sqrt{x+1}+\sqrt{x-1}} \ dx$$ Answers without trigonometric or hyperbolic functions.
Let $F$ denote the function given by this integral defined on a domain with $x > 1$.
Start by noticing the factorization $\displaystyle \frac{2}{\sqrt{x+1}+\sqrt{x-1}} = {\sqrt{x+1}-\sqrt{x-1}}$. Substituting for the denominator of the integrand and simplifying, we have $\displaystyle F = \int 1 - \sqrt{x^2 - 1} \ dx = x - \int \sqrt{x^2 - 1} \ dx$
Next, we recall the integral identity $\displaystyle \int f(x) \ dx = xf(x) - \int xf'(x) dx$ which holds for $f \in C^1\!$. Applying this to $f = \sqrt{x^2-1}$ we get $\displaystyle F = x - \int \sqrt{x^2 - 1} \ dx = x - x\sqrt{x^2-1} + \int \frac{x^2}{\sqrt{x^2-1}} \ dx$

Write $\displaystyle \frac{x^2}{\sqrt{x^2-1}} = \sqrt{x^2-1} + \frac{1}{\sqrt{x^2-1}}$ so we can solve for the integral of $f$:

$\displaystyle 2 \int \sqrt{x^2-1}\,dx = x \sqrt{x^2-1} - \int \frac{1}{\sqrt{x^2-1}}\,dx$

Finally, we can solve this last integral by trig substitution with $x = \sec (s)$ and noting that $\int \sec s \ ds = \log( \tan s + \sec s)$.
Adding up this mess, we get $\displaystyle F = x - \frac{1}{2} \left(x \sqrt{x^2-1} - \log{\large(}\sqrt{x^2-1} + x{\large)} \right)$.

Sorry about the formatting. For some reason, only inline math mode will compile properly.

Anyways, I remembered the following integral which kept me busy for a few hours (or more) in grad school. Enjoy.

Compute
$\int_0^{2\pi} \log \Gamma (\frac{x}{2 \pi}) \exp(\cos x) \sin(x + \sin x) \ dx$ where $\Gamma$ is the usual gamma function.

All times are GMT -8. The time now is 02:56 AM.


text:
# Line spacing with mixed font sizes in itemize environment

I'd like to mix in footnote-sized text with normal-sized text in beamer. The problem is that when I use \footnotesize the line spacing still corresponds to \normalsize spacing. Here's a minimal sample document demonstrating the problem:

\documentclass{beamer}
\begin{document}
\begin{frame}{Line spacing for footnote-sized text in itemize environment}
Problem: footnote-sized text still has normal-sized line spacing.
\begin{itemize}
\item The problem is below:\\
\footnotesize This is some footnote-sized text. When it wraps around it still has normal-sized line spacing.
\normalsize
\item Some other items for comparison
\item I would like the footnote-sized text to have single spacing relative to the footnote size.
\end{itemize}
\end{frame}
\end{document}

I was hoping to use the footnote-sized text to write parenthetical comments in a way which visually distinguishes them from the main items in the list.

-

## migrated from stackoverflow.comApr 26 '12 at 23:03

This question came from our site for professional and enthusiast programmers.

Welcome to TeX.sx! Your question was migrated here from Stack Overflow. Please register on this site, too, and make sure that both accounts are associated with each other (by using the same OpenID), otherwise you won't be able to comment on or accept answers or edit your question. –  Werner Apr 26 '12 at 23:09
Try putting \par immediately before \normalsize. See this answer for details. –  Ian Thompson Apr 26 '12 at 23:14

\documentclass{beamer}
\begin{document}
\begin{frame}{Line spacing for footnote-sized text in itemize environment}
Problem: footnote-sized text still has normal-sized line spacing.
\begin{itemize}
\item The problem is below:

{\footnotesize This is some footnote-sized text.
When it wraps around it still has normal-sized line spacing.\par}

\item Some other items for comparison
\item I would like the footnote-sized text to have single spacing relative to the footnote size.
\end{itemize}

{\footnotesize This is some footnote-sized text.
When it wraps around it still has normal-sized line spacing.\par}

\end{frame}
\end{document}

When a font size changing command is issued, a \par (or an empty line) before the end of the text should be used.

-

I fiddled a little big and am not exactly happy but can share what I found. I tried the following:

\item The problem is below:\\
\begin{spacing}{0}
\footnotesize This is some footnote-sized text. When it wraps around it
still has normal-sized line spacing.
\end{spacing}

This worked okay. I assumed I'd just need to fiddle with the number argument of the spacing environment, but larger than 0 just seemed to leave the compact line spacing of the footnotesize paragraph the same while increasing the gap between that line and the next item. It looked "off." I also tried this:

\item The problem is below:\\
\footnotesize This is some footnote-sized text. When it wraps around it
still has normal-sized line spacing.

This reduces the spacing and the number can be fiddled with but it doesn't seem to want to take consecutive arguments; once linespread is set, it didn't want to set it back to 1 for the subsequent items.

Sorry to not be of more help. I did my digging primarily here:

Update 1/6 If italics is okay, this looked favorable to me?

\item The problem is below:\\
\footnotesize \begin{quote}This is some footnote-sized text.
When it wraps around it still has normal-sized line spacing.\end{quote}

For some reason the quote environment seems to undo the spacing issue and the lines look compressed to me. What do you think? I tried to find out how to define the quote environment in order to remove the italics or indent but was unsuccessful so far.

-
Hi Hendy, thanks for this. Using setspace.sty (for the spacing environment) seems to make all of my footnotes invisible in beamer. Should I start this as a separate question? –  Henry Jan 5 '11 at 10:15
That's odd... not sure what to say about that. See my answer for one more idea? –  Hendy Jan 6 '11 at 16:05

You can undo the italics in the quote environment by adding \em inside it. (Regular type is the italic's emphasis.) To change the quote environment to use regular type globally, use \setbeamerfont{quote}{shape=\sf}.

-


text:
# zbMATH — the first resource for mathematics

$$L^{\infty}$$ error estimates and superconvergence results for a collocation-$$H^{-1}$$-Galerkin method for elliptic equations. (English) Zbl 0584.65073
The paper deals with a mixed approximation scheme based on a collocation method and an $$H^{-1}$$ Galerkin procedure for solving general elliptic equations of the form $$\nabla \cdot (a\nabla u)+b\cdot \nabla u+cu=f$$ in $$\Omega$$, $$u=0$$ on $$\partial \Omega$$. The author establishes $$L^{\infty}$$-error estimates and proves that at particular points in the domain, the rate of convergence is higher than the global rate, i.e. some superconvergence phenomena occurs. Some numerical results are presented.
Reviewer: J.Haslinger

##### MSC:
 65N35 Spectral, collocation and related methods for boundary value problems involving PDEs 65N30 Finite element, Rayleigh-Ritz and Galerkin methods for boundary value problems involving PDEs 65N15 Error bounds for boundary value problems involving PDEs 35J25 Boundary value problems for second-order elliptic equations
Full Text:


text:
• # Simple but Effective Collisions Part 1: Radial Collision Handling

Math and Physics

## The Concept

As mentioned, this system allows you to handle collisions with objects which can be enclosed by bounding cylinders. These cylinders can be defined by two numbers, or rather, a vector and a float - one containing a location, and the other containing a radius. Let's imagine we have a simple scenario with a triangle of trees placed on a simple flat plane:
The player can move the camera around in a first person view, and is able to walk among the trees. We will use the 'radial' collision detection method to prevent the camera from passing through the trees. When the player attempts to walk through a tree, the algorithm will slide them smoothly around it, as seen in pretty much all current first person games. This process involves two steps, which are repeated every frame:
1. If the camera has entered a tree's bounding cylinder, register a collision.
2. Calculate a new location for the camera which puts it back outside the bounding cylinder.
So, how do we do this in practice?

## The Implementation

We'll start out by defining the location of each tree using an array of 2D vectors. We only need 2D vectors because the trees are not going to affect our motion in the up-down direction. If you wanted to use this algorithm for objects on 3D terrain rather than a 2D plane, it should work just as well with the same vector format because it's not possible for a walking person to collide with a tree from above or below, only from the side. Because of this, we only need to worry about the positions of the trees in the x-z plane, whether we have a flat terrain or not. As mentioned, we will use DirectX syntax and functions from here on in:  //The array contains 3 vectors - one for each tree: D3DXVECTOR2 v_treePositions[3] = {D3DXVECTOR2(0, 1), D3DXVECTOR2(1, 0), D3DXVECTOR2(-1, 0)};  We will also define a variable to describe the radius of each tree trunk. We will assume they're all the same radius, but if you wanted to personalise each one, you could define an array of numbers, or add a third component to the vectors describing tree locations.  //Each tree is 1m in radius: float f_treeRadius = 1.0f;  ...and while we're at it, we will define a vector to hold the location of our first-person camera as it changes from frame to frame:  //This is a global vector holding the camera's location. It starts at the centre of the world (0, 0, 0): D3DXVECTOR3 v_camPos(0, 0, 0);  OK, now each tree has an entry in our program, and we know where our player's camera is from frame to frame. All that remains is to put this information to use! We will define a function which will be called every frame. This function will loop through the tree locations, comparing each one to the location of the camera for that frame. If it registers a collision (i.e. if the camera's location is closer to one of the trees than is permitted by the bounding radius), then we will push the camera back outside the bounding radius. The effect is cumulative, meaning that each tree makes its own contribution to the overall position change. This way, the algorithm also works for closely bound groups of trees where the camera could be in contact with more than one at once. Let's go through the function line by line:  D3DXVECTOR2 handle_Collisions_Radial() { D3DXVECTOR2 v_deltaPos = D3DXVECTOR2(0, 0); for (int i = 0; i < 3, i++) { D3DXVECTOR2 v_vectorToCentre(v_camPos.x - v_treePositions.x, v_camPos.z - v_treePositions.z); if (D3DXVec2Length(&v_vectorToCentre) < f_treeRadius) { v_deltaPos += ((f_treeRadius / D3DXVec2Length(&v_vectorToCentre)) - 1) * v_vectorToCentre; } } return v_deltaPos; }  First, we notice that the function returns a D3DXVECTOR2. This corresponds to the amount by which the camera needs to move to put it outside of the bounding cylinders of all the trees. It returns a D3DXVECTOR2 because the trees only affect the location of the camera in the x-z plane. Second, we create an inline vector called v_deltaPos which will be filled with the overall change in position required to put the camera outside the bounding radii of all the trees. This is the vector which the function will eventually return. Next, we enter a loop which will iterate through all the trees. Since we have three trees, we will have three iterations. Third, we create a temporary vector called v_vectorToCentre. This contains a directional vector between the location of the camera and the location of tree in the x-z plane. Fourth, we compare the length of this vector to the bounding radius. Note that the +y direction is up. Fifth and finally, if the camera is indeed within the bounding radius of the current tree, we will add on just enough of v_vectorToCentre to put it back outside the tree trunk. Once the loop is complete, we return the resulting position change. The process is fairly intuitive until the final step. This diagram explains the maths behind it a little clearer:
As you can see, the location of the camera has moved within the bounding radius of tree . In order to put the camera back outside of the bounding radius while simultaneously making it slide smoothly around the edge, we need to add a certain amount of the vector v_vectorToCentre to the camera's position until it is once again outside the bounding radius, such that: In this equation, |v_vectorToCentre| denotes the length of the vector between the camera position and the tree's centre. It is calculated using the D3DXVec2Length() function. We are looking for the value of k - that is, the linear multiple of v_vectorToCentre which we should add to our camera position in order to put it back outside the bounding radius. In order to find the value of k, we can first divide through by |v_vectorToCentre|. This gives us: This can easily be rearranged to get: We have now calculated our value of k. All that remains is to add this linear multiple of v_vectorToCentre to our overall change in location (v_deltaPos). If you review the fifth step in the function, you can see that this is exactly what is happening. If called every frame, this function will calculate a collision contribution for each tree in the loop. Once this contribution is acquired, it should be added to the camera position vector as follows:  //Call the function to assess collisions: D3DXVECTOR2 v_camPosCollisionDelta = handle_Collisions_Radial(); //Since this function returns a D3DXVECTOR2 and our camera position is a D3DXVECTOR3, we //will need to create a dummy D3DXVECTOR3 with no y-component to add to v_camPos: D3DXVECTOR3 v_camPosCollisionDelta3D(v_camPosCollisionDelta.x, 0, v_camPosCollisionDelta.y); //Finally, we will move the camera according to our collison delta vector: v_camPos += v_camPosCollisionDelta3D; //Now that we have our new camera position, we can continue to render the frame as normal...

# Conclusion

That's it! You may be surprised at how well this simple method works. The effect is quite satisfying, and it certainly added an extra professional depth to my exploits in the first-person camera world. Of course, this method could be used to handle collisions involving objects other than a first-person camera, and subjects other than trees. This is also a bare-bones version of the algorithm, and it could easily be adapted to handle collisions in 2D games as well as collisions between circles and spheres and their environment (throughout, we assumed that the first person camera was a simple point rather than an object with dimensions). I intend to write another article describing another collision method which can be used to simulate walls and interiors/exteriors of buildings and things like hedges and maze walls. Until then, thanks for reading!

# Article Update Log

19 May 2013: Initial release

Report Article

## User Feedback

You need to be a member in order to leave a review

## Create an account

Register a new account


text:
# Epstein frame

Epstein frame

An Epstein frame or Epstein square is a standardised measurement device for measuring the magnetic properties of soft magnetic materials, especially used for testing of electrical steels.

The International Standard for the measurement configuration and conditions are defined by the standard IEC 60404-2:2008 Magnetic materials - Part 2: Methods of measurement of the magnetic properties of electrical steel sheet and strip by means of an Epstein frame[1] published by International Electrotechnical Commission.

An Epstein frame comprises a primary and a secondary winding. The sample under test should be prepared as a set of a number of strips (always a multiple of four) cut from electrical steel sheet or ribbon. Each layer of the sample is double-lapped in corners and weighted down with a force of 1 N[clarification needed] (see photo).

The power losses are measured by means of a wattmeter method in which the primary current and secondary voltage are used. During the measurement, the Epstein frame behaves as an unloaded transformer.

Power loss, Pc, is calculated as:

${\displaystyle P_{c}={\frac {N_{1}}{N_{2}}}\cdot P_{m}-{\frac {\left(1,111\cdot |{\bar {U_{2}}}|\right)^{2}}{R_{i}}}}$

where:
${\displaystyle N_{1}~}$ is the number of turns of primary winding
${\displaystyle N_{2}~}$ is the number of turns of secondary winding
${\displaystyle P_{m}~}$ is the reading of the wattmeter in watts
${\displaystyle R_{i}~}$ is the total resistance of the instruments in the secondary circuit in ohms and
${\displaystyle |{\bar {U_{2}}}|}$ is the average secondary voltage in volts.

Specific power loss, Ps, is calculated as:

${\displaystyle P_{s}={\frac {P_{c}\cdot 4\cdot l}{m\cdot l_{m}}}}$

where:
${\displaystyle l~}$ is the length of the sample in metres
${\displaystyle l_{m}~}$ is the average magnetic path length = 0.94 (constant value)
${\displaystyle m~}$ is the mass of the sample in kilograms

If all conditions are as defined in the standard, the standard deviation of the reproducibility of the values is not greater than 1.5% up to 1.5 T[clarification needed] for non-oriented electrical steel and up to 1.7 T for grain-oriented electrical steel.


text:
# Tower of Hanoi

Rings:

Steps: 0

## About the Tower of Hanoi

The Tower of Hanoi is a game created by Édouard Lucas, a French mathematician, in 1883. The object of the game is to move the stack of $$n$$ disks to another rod, in their original order. Only one disk may be moved at a time, and a disk may never be placed on top of a smaller disk. The minimum number of moves required in any game is $$2^n - 1$$.

## Using the Applet

Choose the number of disks you would like to start with from the menu. Click on a stack to select its top disk, then click the stack you would like to move that disk to. You cannot place a disk onto a stack with a smaller disk. To unselect a disk, click on its stack again.

The text field next to the game displays a binary string that tracks the movements of each disk. The $$i$$ bit in the string refers to the $$i$$ disk, where the largest disk has index 0. When a disk is moved, its state is changed. Note that the sequence of strings does not encode where the disk was placed nor the current state of the game at any time; only the identity of each disk moved is recorded. One of the optimum solutions to the puzzle will yield the Gray code.


text:
# How do you multiply 6h+5(h-5)?

Jun 30, 2018

$11 h - 25$

#### Explanation:

$6 h + 5 \left(h - 5\right)$

First, use the distributive property (shown below) to simplify $5 \left(h - 5\right)$:

Following this image, we know that:
$\textcolor{b l u e}{5 \left(h - 5\right) = \left(5 \cdot h\right) + \left(5 \cdot - 5\right) = 5 h - 25}$

Put it back into the expression:
$6 h + 5 h - 25$

Add $\textcolor{b l u e}{6 h}$ and $\textcolor{b l u e}{5 h}$:
$11 h - 25$

Hope this helps!


text:
# More Friction

1. Nov 22, 2005

### Izekid

Here it comes, I have a frictiongrade who is 0,45 A car who is driving in 18m/s and mass of the car 1200kg

For the Frictionpower it is : 0,45 * (mg) = 5,5 kN

And how do I calculate the meters the car travels until it stops?

2. Nov 22, 2005

### skeeter

initial kinetic energy - work done by friction = 0

(1/2)m(vo)^2 - f(delta x) = 0

solve for delta x

3. Nov 23, 2005

### andrevdh

Use the work-kinetic energy theorem. That is the work done by all of the forces acting on an object (or the resultant of all of these) is equal to it's change in kinetic energy:
$$W_{F_R}=\Delta K$$


text:
Browse Questions
Home  >>  AIMS

Which of the following molecular formulae belongs to alkyne series?

$\begin{array}{1 1}(A)\;C_7H_{14}&(B)\;C_{10}H_{22}\\(C)\;C_9H_{16}&(D)\;C_{16}H_{12}\end{array}$

$C_9H_{16}$
Hence (C) is the correct answer.


text:
# Important Programming Concepts (Even on Embedded Systems) Part II: Immutability

September 14, 20141 comment

Other articles in this series:

This article will discuss immutability, and some of its variations in the topic of functional programming.

There are a whole series of benefits to using program variables that… well, that aren’t actually variable, but instead are immutable. The impact of immutability on software design can differ quite a bit, depending on the type of computing platform, and especially whether dynamic memory or concurrent memory access are involved.

When I took my first computer science course in college, the instructors extolled the virtues of immutable data. Here’s the message of the Functional Weenies Wizards: Functions that alter data are flawed and evil, and must be marked (e.g. the exclamation point at the end of set! in Scheme), like the Mark of Cain or the Mark of the Beast; functions that do not alter data are pure and good. I thought they were kind of extreme. (The FW is the guy in the corner at the party who walks up to you and says “Hey, I’ll bet you didn’t realize you’ve been using monads all along!” And when you hear the word “monad”, that’s your cue to run away. Or at least play the Monad Drinking Game and down another shot of Scotch.)

But in the intervening years, I’ve come around a bit to their view, and have become a Functional Wannabe. Yeah, mutability has its place, but there are some real advantages to striving towards immutability, and again, it really depends on the type of programming environment you’re in: if the cost of leaving the input data untouched and returning new output data is high, you’re often stuck with mutable data; otherwise, you really should make data immutable wherever possible.

This assumes you have a choice. Let’s forget for a moment about how you can make your programs use immutable data, and first just focus on the difference.

## The Immutable and the Mutable

Consider two pieces of data, the first immutable and the second mutable.

1. The text of the Declaration of Independence
2. The real-time contents of a stock portfolio, namely an unordered list of stock symbols, each with the number of shares owned.

Here’s what you can do with the Declaration of Independence, and what you can’t do with the stock portfolio:

• Use copies of it rather than the original. If all we care about is the text, we don’t care whether it’s in one area of memory or another. If there are seven areas of our program that need access to the text of the Declaration of Independence, we don’t care if all seven refer to one chunk of memory, or three of them refer to one chunk and the other four another chunk, or if each has its own copy. With the stock portfolio, if we make a copy, the copy is going to be out of date and will not be updated if there are any changes. Any piece of software that wants access to the stock portfolio must get it from its source.

• Use the original rather than copies. If no one can change the Declaration of Independence, it’s safe to give out references to it all day. (“Hey bud, wanna see the original copy of the Declaration of Independence? Just go to this bus station locker and look inside; here’s the combination....”) Whereas with the stock portfolio, if we want to give read-only access to other software programs, we have to be paranoid and create a copy, so that we can prevent malicious programs from altering the original.

• Share access to it without any worry about the problems of concurrency. We can let a thousand different software routines access the text of the Declaration of Independence concurrently. It doesn’t matter. None of them are going to be able to change it, so there’s no risk of data corruption. The stock portfolio is another story. We can let only one software task at a time modify it, and while that modification is taking place, we can’t allow other software tasks to read it.

• Simplify expressions at compile time. In compiled software, the compiler can precompute the number of words in the Declaration of Independence; we don’t have to execute operations at runtime to compute the number of words. But we can’t count the number of stocks in the portfolio at compile time, because it can change.

• Cache the result of expressions computed at run-time. We could compute a million different statistics about the Declaration of Independence: the number of times the letter pairs ou, ha, or ty occur; the list of words in alphabetical order that score less than 12 points in Scrabble; the complete Markov chain transition matrix of two-letter sequences in the Declaration; the MD5 and SHA1 hashes of the full text; the ROT13 translation of the all the words with at least 10 letters… and so on. We can’t possibly compute them all in advance at compile-time, so it’s not worth doing. But if we do happen to compute some summary statistic once, we can save the results for later, because the Declaration of Independence will never change. We can’t do that with the stock portfolio though.

• Store it in nonvolatile memory. In embedded microcontrollers, the amount of available program memory might be 5 or 10 times larger than the amount of available RAM. This means that RAM is more precious than program memory. We can store the text of the Declaration of Independence in program memory, leaving that precious RAM for something else. The stock portfolio is mutable, so we can’t store it in program memory.

• Use it as a key in an associative array. Associative arrays are ubiquitous in software. They’re just sets of key-value pairs. Use the key to index a value. There are all sorts of data structures to do so efficiently. The only catch is that the keys have to be immutable (or at least comparable and/or hashable in an immutable way), otherwise we may not be able to find the value again.

## Programming Language Support for Immutability

In some programming languages there are keywords that specify a degree of immutability. These are extremely useful, but you should be careful that you know exactly what these keywords mean if you do use them. C/C++ have const and Java has final, but these have different meanings, and neither one marks an object as truly immutable.

In C and C++, the const keyword has several usages. When applied to a program variable or function argument, it means that the code given access to that data is not allowed to change it, within the appropriate scope:

int doofus123(const data_t *mydata, data_t *mydata2)
{
const data_t *mydata3 = mydata2;
...
}

Here the variables mydata and mydata3 are marked as pointers to const data. The code within doofus123() is not allowed to change the data referenced through these pointers. But that doesn’t mean that this data is guaranteed not to change in other ways. The doofus123() function could change the data through the pointer mydata2, which points to the same location as mydata3. And we can’t assume that mydata and mydata2 point to different locations (if they do refer to the same data, this is called pointer aliasing), so there’s no guarantee that the data referenced by mydata won’t change.

In C++, if an object’s method is marked with const, it means that method cannot modify any of the non-mutable fields of the object, and that method can be called through a const * to the object, whereas non-const methods cannot be called via a const *:

class Dingus
{
int data;

public:
void setData(int x);
int getData() const;
};

Dingus dingus;
const Dingus *pdingus = &dingus;
int n = pdingus->getData();   // OK
pdingus->setData(33);         // illegal: pdingus is a const *

You should learn to use const! It’s the key to writing software that promises not to change data. The compiler will produce an error to stop you from inadvertently violating this promise, at least if you don’t circumvent the compiler and cast away the const qualifier. In an embedded system, it also allows you to mark initializer data in such a way that it can be stored in nonvolatile program memory rather than in generally scarcer volatile data memory:

const char author[] = "James Q McGillicuddy";

Neither the const * concept in C or the const & concept in C++ (data that the compiler will not allow you to change through const pointers and references) has an equivalent in Java; the Java final keyword applies only to the direct contents of a variable itself, and it just means that the Java compiler will not allow you to assign a final field more than once.

public class Doofus
{
private int x;

public void setX(int newX) { this.x = newX; }
}

public class Dingus
{
final private Doofus doofus;

public Dingus()
{
this.doofus = new Doofus();
}
public void tweak()
{
this.doofus = new Doofus();  // illegal: doofus is final
this.doofus.setX(23);        // OK, we can call any public method
}
}

Neither language provides a way to mark data as truly immutable, however; we can only hide that data behind an API that prevents cooperating software from modifying the data.

The poster child in Java of immutability, or lack thereof, is the mutable java.util.Date class. If you want to use Date objects to record the date and time of your birthday or anniversary or upcoming colonoscopy, and you allow someone else access to those objects, you are giving them a blank check to change the fields of the Date object. To be safe, you need to make a defensive copy. That’s extra work that would be unnecessary if Date had immutable fields.

In fact, there’s a whole school of thought that everything should be immutable; the functional programming languages strongly encourage the use of immutable objects. When you want to change some object X, instead of altering its data, you create a whole new object X in its place. This sounds like a lot of extra work — but anytime you think that, go look at the advantages I cited earlier with immutable data.

## An Example: The Immutable Toaster

The embedded guys who deal with plain old integer data types in C are probably looking at this and scratching their heads. Objects don’t enter the picture at all, and “new object” seems to allude to dynamic memory allocation, which a lot of us in the embedded world avoid like the plague. So here’s an example; let’s write a program that controls a toaster. It has an update function that’s called 10 times a second. Oh, and we need to modulate the heating element on and off, keeping it on at most half the time, not all the time, otherwise it will be too hot and have a tendency to burn the toast.

enum ToasterState { TS_INACTIVE, TS_TOASTING, TS_READY };

typedef struct TAG_ToasterInfo {
ToasterState state;
int16_t timer;
bool heating_element_on;
} ToasterInfo;

void toaster_update(ToasterInfo *pti, bool buttonpressed, bool dooropen, int16_t dial)
{
if (pti->state == TS_INACTIVE)
{
if (buttonpressed)
{
pti->state = TS_TOASTING;
pti->timer = 0;
pti->heating_element_on = true;
}
}

if (pti->state == TS_TOASTING)
{
if (pti->timer >= dial)
{
pti->heating_element_on = false;
}
else
{
++pti->timer;
if (!pti->heating_element_on)
pti->heating_element_on = true;
if (pti->heating_element_on)
pti->heating_element_on = false;
}
}

{
if (dooropen)
{
pti->state = TS_INACTIVE;
}
}
}

Very simple, you turn the dial to set the toasting time, press the button, the toaster turns on and starts a timer, toggling the heating element on and off, and when the timer exceeds the dial setting then the toaster turns off and waits for you to open the door before it will start again. Not much to it!

Can you find the bugs here? They’re subtle.

One problem is that we read and write pti->state all over the place in this function, and we muddle its new and old value. By “muddle”, I mean that we update pti->state at one point in the function, and then we read its value later, when we probably meant to refer to its original value at the start of the function. In this case, the effect is not that critical; it just means the toaster finishes what it’s doing a split second earlier, depending on the order of statements within the function; if we rearrange them, we might get state transitions that occur one iteration later.

The bigger problem is in the heating element toggle code:

     if (!pti->heating_element_on)
pti->heating_element_on = true;
if (pti->heating_element_on)
pti->heating_element_on = false;

This turns heating_element_on to true, but then immediately afterward, turns it to false. Again, we are muddling its new and old value.

You might think this is an obvious bug, and no sane person would write code that works that way. Why not just toggle heating_element_on:

pti->heating_element_on = !pti->heating_element_on;

Well, we could, in this case, and it would be correct. But that’s only because this is a really easy kind of behavior to design. What if we needed to do something more complicated, like run the heating element at 60% duty cycle, or if we needed to control the heating element based on a thermostat?

The problem with changing state variables in place, is that we have plenty of opportunities to write programs with errors, because the same name in a program now refers to different values, which are sensitive to the order in which we do things. Let’s say there’s some other logic in the program that needs to do something at the precise moment when the heating element is on and something else occurs. If this other logic precedes the toggle statement, then pti->heating_element_on refers to the previous heating element state, whereas if the other logic follows the toggle statement, then pti->heating_element_on refers to the next heating element state. And that’s a simple case, because there’s only one place in which the heating_element_on member is modified. In general it’s not that simple; there may be multiple places in the code where state can be modified, and if it occurs in if statements, then sometimes the state is modified and sometimes it is not. When we see pti->state or pti->heating_element_on, we can’t tell if these refer to the value of state variables as they were at the beginning of our update function, or a changed value at the end, or some intermediate transient value in the middle of calculations.

(Things get even uglier if there are complex data structures, like a linked list, where a pointer itself is modified:

plist->element.heating_element_on = true;       // Line A
plist = plist->next;
plist->element.heating_element_on = false;      // Line B

Here we have two lines, A and B, which both change mutable state, and the lines look the same, but they refer to two separate pieces of data.)

A better design for our toaster software uses separate data for input and output state:

void toaster_update(const ToasterInfo *pti, ToasterInfo *ptinext,
bool buttonpressed, bool dooropen, int16_t dial)
{
/* default: do what we were doing */
ptinext->state = pti->state;
ptinext->timer = pti->timer;
/* for safety: default the heating_element_on to false */
ptinext->heating_element_on = false;

if (pti->state == TS_INACTIVE)
{
if (buttonpressed)
{
ptinext->state = TS_TOASTING;
ptinext->timer = 0;
ptinext->heating_element_on = true;
}
}

if (pti->state == TS_TOASTING)
{
if (pti->timer >= dial)
{
ptinext->heating_element_on = false;
}
else
{
ptinext->timer = pti->timer + 1;
if (!pti->heating_element_on)
ptinext->heating_element_on = true;
if (pti->heating_element_on)
ptinext->heating_element_on = false;
}
}

{
if (dooropen)
{
ptinext->state = TS_INACTIVE;
}
}
}

To use this version of the update function properly, we need to avoid pointer aliasing, so that inside the function, we have complete freedom to change the contents of the next state, while still being able to assume that the contents of the existing state stays unchanged. So we can’t do this:

ToasterInfo tinfo;
tinfo.state = TS_INACTIVE;
while (true)
{
/* get inputs here */
toaster_update(&tinfo, &tinfo, buttonpressed, dooropen, dial);
/* set output here based on tinfo.heating_element_on */

wait_100_msec();
}

But we could call the update function using a pair of ToasterInfo states:

ToasterInfo tinfo[2];
ToasterInfo *ptiprev = &tinfo[0];
ToasterInfo *pti     = &tinfo[1];
ptiprev->state = TS_INACTIVE;
while (true)
{
/* get inputs here */
toaster_update(ptiprev, pti, buttonpressed, dooropen, dial);
/* set output here based on pti->heating_element_on */

/* swap pointers to state variables */
ToasterInfo *ptmp = pti;
pti = ptiprev;
ptiprev = ptmp;

wait_100_msec();
}

Now when we are looking at the toaster_update() function, we can clearly distinguish the old value of the toaster state from the new value of the toaster state. And because we used const, we can have the compiler help catch our errors if we try to assign to any of the pti members instead of ptinext. Moreover, if we have the program stopped in a debugger, we can see both the previous state and the new state in their entirety, without having to do any clever sleuthing to infer this information.

This creates a little bit of extra work, but one of the key lessons of software design is that you have to make tradeoffs, and often you end up giving up a little bit of performance efficiency to gain improvements in code clarity, modularity, or maintenance.

This kind of approach (e.g. new_state = f(old_state, other_data)), where we decouple the mutable state update from the computation of new state, is something I call “pseudo-immutability”. It’s not purely immutable, since we are storing a mutable state variable somewhere, but the data looks immutable from the standpoint of computing the new state, and the software design tends to be cleaner as a result.

## A Tour of Immutable Data

One key idea when considering immutability is the value type. Values are just ways of interpreting digital bits, and the important point is what they are, not where they are stored.

Let’s say you have a software program that works with a complex number e.g. 4.02 + 7.163j; this is just a way of interpreting two 64-bit IEEE-754 floating-point values, one for the real part, and one for the imaginary part. This number exists independently of where it is stored. In fact, if we want to be purely functional, there is no storage: there are only inputs and outputs.

When we do have state containing a complex number, it happens to be just a transient place to store the 128 bits. The complex numbers themselves are nomads (“Hey, I’ll bet you didn’t realize you’ve been using monads all along!” “No, I said nomads, not monads, you creep!”) traveling from place to place, stopping inside a state variable only for a microsecond. One problem with the average approaches of programming using mutable state, is that this idea of transiency disappears. Instead, we get focused on some variable z containing a complex number, that lasts for a long time, and we poke and prod at it, and the name z might refer to an erratically changing series of values.

The alternative to a value type is a reference type, where a critical aspect of the data is where it is stored. If I have persistent mutable state in a software program, like a color palette for a window manager, the data lives somewhere, and unless I want to deal with potentially out-of-date copies, when I work with that color palette, I pass around a pointer to that data. The value of the color palette is a particular fixed choice of colors. But the window manager’s reference to the color palette is a unique container for those values; it lives in one place, and its contents can change. (In C, the term lvalue refers to a storage location, which, if it does not have a const qualifier, can be used on the left side of an assignment statement; the term rvalue refers to an expression which can be used on the right side of an assignment statement.)

### Microarchitecture: Low level immutable data

Perhaps when you program, you picture RAM as kind of like one of those plastic storage containers with lots of little compartments. Locations 0x3186-0x319c contain your name; locations 0x319d-0x31f2 contain your address, and so on. At a low level, this is because RAM retains state by default: as long as we keep it happy, some pattern of bits at a particular location will stay the way they are. With static RAM (SRAM), we just provide a well-conditioned source of power; with dynamic RAM (DRAM), we have to activate refresh circuitry, which used to be part of the overall system design in a computer, but is now an intrinsic part of the memory module itself. In either case, the RAM requires an explicit write step to change its contents. On the other hand, if you look at the registers and data paths inside a processor, they are very functional in nature. Here’s the block diagram of the Microchip PIC16F84A microcontroller, from the datasheet:

The 16F84A has a register file for its data memory. The difference between a register file and SRAM is kind of a gray area that depends on the architecture and implementation. The classical definition of a register file is a small group of registers (one for each address in the register file) with separate input and output data ports; we can think of a register as a tandem group of D flip-flops, which are clocked data storage devices: during a particular clock cycle, the output of a flip-flop is constant, but at the rising edge of the clock signal, the D flip-flop’s input is captured by the flip-flop and propagated to the output during the next cycle. Memory in a D flip-flop is transient and only lasts for 1 clock cycle. If you want the data to persist in value, you have to feed the output back to the input. (Those of you who are more familiar with DSP than digital design can think of a register or a D flip-flop as a unit delay z-1.)

So a register file can be thought of as a bunch of D flip-flops and multiplexers: If you are writing new information to the register file to a particular location, the multiplexer for that register gets its data from the register file’s input port, whereas if you’re not writing to the register file, the multiplexer maintains each register’s value by getting its data from the previous value of the register.

The special-purpose registers of the 16F84A, like the program counter or the W register, are similar, but they have manually-designed data paths:

The W register is an 8-bit register which is always a function of its previous value and some other data. The exact function used depends on the arithmetic-logic unit (ALU). And this pattern is essentially the same as the “pseudo-immutable” approach I have been talking about:

 new_state = f(old_state, some_other_data);

This keeps new and old state conceptually separate, as opposed to the average old mutable approach that causes us headaches:

 // "state" refers to old values
if (some_condition)
{
state.some_member = state.some_member + something_else;
}
// now "state" refers to new values
if (other_condition)
{
state.other_member = state.some_member - another_thing;
}
// now "state" has changed again

So it’s not just the Functional Wizards who are using immutable data; at a low level, many microprocessors and microcontrollers utilize the same approach of separating input and output data.

### Immutable Data Structures

Back to the high level of programming: there’s a whole subject of immutable or persistent data structures. The idea is that you don’t change the contents of data; instead, you create a new set of data.

This may seem wasteful and impractical, and in some cases it is. For example, if you have an array X of 1 million integers, and you want to change the 30257th element of X, it takes a lot of work just to create a new array with a new value for element #30257 and copies of all the other elements of X. Then if you don’t need the old array any more, you throw it away. Wasteful… except that general-purpose computers are fast and have lots of memory, and the high-level languages recycle unused memory through garbage collection. In Python, the numpy library uses optimized C code underneath to manipulate arrays, and although mutable data access is supported, there are many numpy methods that are functional in nature and create new copies of data rather than changing state in-place. The advantages usually outweigh the costs; we can design programs with fewer errors, and by separating input and output we can optimize the computation, by using parallel computations or a GPU to speed up functional programming, whereas mutating state prevents many of these optimizations.

But arrays are only one type of data structure.

For most of the usual mutable data structures, there are immutable equivalents. The classic example is the list. A mutable list might rely on a linked list implementation. Immutable lists also use linked lists, but they have that nice property, like any other immutable data structure, that the nodes of an immutable list won’t change. Here’s an example in Java:

 public class ConsList<T>
{
private final T car;
private final ConsList<T> cdr;
public ConsList(T car, ConsList<T> cdr)
{
this.car = car;
this.cdr = cdr;
}
public T head() { return this.car; }
public ConsList<T> tail() { return this.cdr; }
}

Really simple: a list consists of one item (the head() of the list) and a reference to the rest of the list (the tail()). The names cons and car and cdr are weird names that date back to old LISP implementations, and perhaps it would have been better to get rid of cons, and use first and rest:

 public class ImmutableList<T>
{
private final T first;
private final ImmutableList<T> rest;
public ImmutableList(T first, ImmutableList<T> rest)
{
this.first = first;
this.rest  = rest;
}
public T head() { return this.first; }
public ImmutableList<T> tail() { return this.rest; }
}

If we wanted to store the list (1,2,3,4) we would do it this way:

ImmutableList<Integer> list1 =
new ImmutableList<Integer>(1,
new ImmutableList<Integer>(2,
new ImmutableList<Integer>(3,
new ImmutableList<Integer>(4,
null))));

We can visualize the list nodes as follows. Each contains a pair of references, one to the list head and the second to the remaining elements, with the end of the list denoted by a link to null:

Insertion at the beginning is easy:

ImmutableList<Integer> list2 =
new ImmutableList<Integer>(5, list1);
// this is the list (5,1,2,3,4)

Decapitation (removal of the head) is also easy:

ImmutableList<Integer> list3 =
list1.tail();
// this is the list (2,3,4)

While these three lists are conceptually separate, they reuse many list nodes.

In fact, the reason these are called persistent data structures, is because the immutability of the data means that you can keep around references to older versions, and they will be unaffected by the changes to the data structure, since really there are no changes, only new structures made of nodes which may be newly allocated or may be shared with other immutable data.

Other operations (appending, or insertion/removal of interior nodes) are possible but require extra list node copies and take longer to finish. In Java, this type of list is a bit verbose to manipulate; in functional languages like LISP or Haskell, list node concatenation is a primitive operation that is much shorter to write in code.

Other data structures like queues, maps, trees, are also possible to implement using immutable techniques. There has been a great deal of research on efficient implementations of these structures, including so-called “amortized algorithms” for reducing the worst-case execution time for things like binary tree rebalancing, by spreading the occasional long operation into many short operations that extend into the future. It’s just like a home mortgage: instead of having to come up with $399,000 in one fell swoop, we go to the bank and get a 30-year 4.5% mortgage with a monthly payment of$2001. We can keep our worst-case costs per operation low by spreading them out over time.

The catch with this whole technique — whether you use complex amortized algorithms, or a simple list with cons cells — is that you need to use dynamic memory allocation, and you need to have some way of recycling unused data, namely reference counting or garbage collection. And that means that while it works great for high-level languages, in low-level languages like C or even C++, immutable data structures are difficult to manage. In C++ there are possibilities (C++11 introduced shared_ptr<> to facilitate automatic memory management; before C++11 the Boost libraries can be used for their own shared pointers); you can read about some of them in these articles.

So immutable data structures are probably out of the running in the low-level embedded world, at least until someone invents a technique for automatic memory management, that also satisfies the hard real-time, determinism, and/or memory safety requirements of many applications.

### Loops?

We’ve talked about the pseudo-immutable approach of isolating updates of mutable state from the computation of new state:

new_state = f(old_state, some_other_data);

There are other ways of enhancing immutability in software design. One of them involves control flow. Take, for example, the simple for-loop in C:

for (int i = 0; i < L; ++i)
{
do_something(i);
}

Here we have a mutable loop variable i. C++ uses a similar approach with STL iterators to facilitate iteration over a generic container, by overloading the increment operator ++ and the pointer dereferencing operator *:

std::vector<Something> mydata = ...;
for (std::vector<Something>::iterator it = mydata.begin();
it != mydata.end();
++it)
{
do_something(*it);
}

Again, a mutable loop variable. And there’s nothing really wrong with it, but if you get a warm happy feeling from functional-style programming, you will not get a warm happy feeling from these programming snippets. The STL use of immutable iterators are a little clumsy (::const_iterator vs ::iterator), and there’s always a risk with more complex software that you might accidentally put another ++it somewhere by mistake, and unintentionally increment the loop twice when you really meant just to do it once. The compiler won’t catch your mistake, as it allows you to modify the loop variable in the body of the loop, and it has no way of figuring out what you really wanted.

There are two approaches for removing mutable loop variables. One is recursion. LISP and other functional languages encourage the use of recursion rather than looping. Here’s an example in Java that uses the ImmutableList<T> type I described earlier:

/* Find the sum of a list's elements.
*/
int sum(ImmutableList<Integer> list)
{
if (list == null)
{
return 0;
}
else
{
}
}

Simple! To find the sum of the elements of a list, it’s either 0 for an empty list, or it’s the head of the list plus the sum of the remaining elements. No loop.

So what’s the catch? It’s this line:

return list.head() + sum(list.tail());

We have to use recursion, which means that if a list has N elements, we need to use up N stack frames. This imposes a memory cost for the stack, and means that if we are debugging and put a breakpoint somewhere for a 1000-element list, there might be 1000 stack frames for us to scroll through in the debugger. Yuck. But that’s not fair, because we know we can solve the same problem with an iterative approach that only uses 1 stack frame:

int sum(ImmutableList<Integer> list)
{
int result = 0;
ImmutableList<Integer> p = list;
while (p != null)
{
p = p.tail();
}
return result;
}

The functional languages have a feature called tail recursion in which the compiler or interpreter is smart enough to recognize recursion that doesn’t require additional stack frames; instead of pushing the return address on the stack and then calling the function again, it sets things up to jump back to the beginning of the function, effectively forming an iterative implementation of a recursive solution. Sometimes the compiler or interpreter is smart enough to do everything on their own; other times it has to be helped by writing code that is amenable to tail recursion. Our example isn’t:

return list.head() + sum(list.tail());

The problem is that we have to first make a recursive call to sum(), then add the head of the list. Adding the head of the list is a deferred action that blocks us from performing tail recursion, because there’s still something to do after the recursion completes. Here’s a modification:

/* Find the sum of a list's elements.
*/
int sum(ImmutableList<Integer> list)
{
return sum_helper(list, 0);
}

int sum_helper(ImmutableList<Integer> list, int partial_sum)
{
if (list == null)
{
return partial_sum;
}
else
{
}
}

In this case, we can handle the addition step before we recurse. The recursion is the last thing, therefore we can theoretically just jump back to the beginning of sum_helper() and not have to waste unnecessary stack frames. But it requires the language implementation to perform tail call optimization, and neither C nor Java appear to guarantee that they will do so, although presumably some C compilers will and maybe Java will in a future version.

The approach I like better is the one used in Python, where instead of recursion, we have a for-loop construct where the iteration is managed within the language:

for item in collection:
do_something_with(item)

Here we have to implement collection in a special way so that the language understands how to iterate over its items. (In Python it’s the iterator protocol.) The language will take care of executing the loop body, with the item variable bound to each item of the collection in sequence. If Python supported immutable data, we could make the loop variable item immutable: even though it takes on different items each time around the loop, within each loop iteration it is intended to be constant.

Of course, in our sum example we would still need to use a mutable state variable:

mysum = 0
for item in collection:
mysum += item

The purely functional alternatives are to use tail recursion, or higher order functions like reduce() or fold() to combine the elements of a collection. We just have to come up with a user-defined function that specifies how each element is supposed to be combined with the result.

Does it matter? Well, yes and no. The mutable approach in this case (mysum += item) is really simple, and it’s very clear to the reader what is going on, so there’s no strong reason to use functional techniques. But in more complicated cases, the functional approaches can lead to a more organized approach that helps ensure correctness. Let’s say we wanted to sum the square of all values which are divisible by 3:

mysum = 0
for item in collection:
if item % 3 == 0:
mysum += item*item

That’s not so bad either… but it’s starting to get more complicated. Functional approaches create a pipeline:

mysum = sum(square_map(div3_filter(collection)))

where we’d have to define functions square_map that yields a new list with squared values, and div3_filter yields a new list of only the items which are divisible by 3. We can either write them manually, or use the higher-order functions map() and filter(). In Python we can do this more concisely with list comprehensions or generator expressions:

sum(item*item for item in collection if item % 3 == 0)

If you’re used to comprehensions, this is very clear and concise. If you’re not, it seems weird, and the original combination of for and if statements would be better. They’re essentially equivalent. You say po-TAY-to, I say po-TAH-to. The difference is that in the functional approach, I don’t have to keep track of the intermediate state involved in the accumulation of sums. And if that doesn’t matter, then pick whichever approach you and your colleagues like best. But sometimes, handling all the housekeeping of intermediate state can be rather cumbersome.

Perhaps the simplest example of how intermediate state makes programming more complex is the swap operation. In Python we can use destructuring assignment:

a,b = b,a

Simple and clear! If we could only assign to one value at a time, we’d be stuck with using a temporary variable:

tmp = a
a = b
b = tmp

That’s fairly simple as well, but then I have to stop and think if I did that right, and Doubt enters my mind, and I think that maybe I should have used this logic instead:

tmp = a
b = tmp
a = b

Hmm. (In case you’re wondering, this second attempt is wrong; the previous code snippet is the correct one.)

In any case, there’s a theme here. To help yourself avoid errors, keep things clear and simple, and removing unnecessary mutable state is often a way to keep things clear and simple.

### Shared Memory?

Another reason for considering immutability involves shared memory. With multiple threads or processes, shared memory is a common technique for using state with concurrent processing. But it imposes a cost, namely of using locks or mutexes to ensure exclusive access.

The immutable way of handling shared state is… well, it’s to avoid sharing state at all. No shared state means no locks. And the classic way to get there is to use message passing. Erlang is probably the prime example of a language with built-in message-passing support, and nearly all GUI frameworks use message queues.

The cost of message-passing is maintaining a message queue, but there are some fairly simple algorithms (such as the non-blocking one described by Michael and Scott) which guarantee proper concurrent access to a queue. Also, if we’re going to use message-passing to overcome the need for shared data, we must use value types: we have to send the data directly. We can’t send pointers to data, since that would require the data to be directly addressable by more than one thread or process, and that would require shared data. Any type of reference that can be sent via message passing (like an array index or a lookup key) must be meaningful to both the sending and receiving entity.

Here’s another example: let’s say we had 10 worker threads, and wanted to use each of them to analyze some fragment of data, and boil the analysis down to a number and add it in to a cumulative total. The shared memory approach looks something like this for each worker thread:

partial_sum = compute_partial_sum(data)
try:
lock(shared_sum_lock)
shared_sum += partial_sum
finally:
unlock(shared_sum_lock)

The message passing approach looks like this:

# Worker thread:
partial_sum = compute_partial_sum(data)

sum = 0
if message_tag == 'partial_sum':
sum += contents

With the shared memory approach, we rely on each of the worker threads to operate correctly. If any of them hangs or crashes before unlocking the shared resource, the whole system can deadlock. With the message-passing approach, we don’t have to worry about this. The only way the whole sum can screw up is if the main thread hangs or crashes; otherwise, the worst outcome is that we might not hear from one of the worker threads, and we can deal with this by rewriting the main thread to time out, and either give up or restart a new worker thread on the missing data.

## Pure functional approaches for managing state

There’s no perfect answer for when to use mutable data and when to use immutable data. The extreme case of eliminating all mutable data whenever possible, is something that is used in some functional programming environments. In Haskell, for example, you have to go out of your way to use mutable objects. And that means you’re probably going to have to talk to the guy in the corner about monads. Or talk to one of the other guys, who are less intimidating, about them. There are a lot of good reasons why they are useful, it just takes quite a bit of mental energy to figure out what’s really going on.

return :: a -> Maybe a
return x  = Just x

(>>=)  :: Maybe a -> (a -> Maybe b) -> Maybe b
m >>= g = case m of
Nothing -> Nothing
Just x  -> g x

When you understand them well enough that they’re intuitive for you to use, let me know what you learned.

## Wrap-up

Here’s a summary of what we talked about today:

• Using immutable data simplifies software design, removing problems for shared access and concurrency, and removing the need for defensive copies.
• Some programming languages include support for some aspects of immutability, like const in C/C++ and final in Java, but these generally don’t provide full immutability.
• Pseudo-immutable techniques are a good real-world compromise to reduce the chances of human error, by decoupling the state update from the computation of its new value, rather than changing mutable state in place:

new_state = f(old_state, other_data)

• The registers in low-level microarchitectures use the pseudo-immutable approach, with flip-flops forming the propagation from old state to new state.
• Immutable (aka “persistent”) data structures can support the use of immutability to create versions of lists, maps, queues that have the advantages of immutability.
• Functional techniques of recursion and iteration are alternatives to looping with mutable state variables.
• Message-passing leverages immutability to eliminate the concurrency problems of shared memory.
• You can use monads if you want to turn to pure functional techniques. (“Hey, I’ll bet you didn’t realize you’ve been using monads all along!”) Remember, I said you can use monads. I probably won’t, at least not intentionally.

Whew! That was longer than I had intended. The next article is on volatility and is much shorter.