This Jupyter notebook\n", "To run a cell, press 'shift+enter' while the cell is selected. Before running, make sure you have installed NUMPY, Matplotlib, and H2MM_C.
\n", "\n", "NUMPY is one of the most commonly used python modules in scientific computing, so it might be included already in your install, but if it is not, you can follow the installation guide here: https://numpy.org/install/\n", "\n", "Matplotlib is a very common python module for creating professional looking plots and graphs, and the default in many scientific applications, it may also be included by default in your python installation. If it is not, see the install guide here: https://matplotlib.org/stable/users/installing.html\n", "\n", "\n", "If you are new to Python, and want some tutorials, may I recommend the tutorials on Corey Schafer's YouTube channel. Some of the particular videos I'd recommend are the following: lists tuple and sets, conditionals, loops and iterations, functions.\n", "\n", "Cells can be 'Code' 'Markdown' 'Raw NBConvert' or 'Heading'\n", "- Code cells are the \"meat\" of the jupyter notebook, they are where you put python code you actually want to run\n", "- Markdown cells are for communicating with the reader, they use the github flavor of Markdown, so when they are rendered, (either when you open the notebook, or after editing you can hit shift+enter) the text is readable, and has various formatting niceties like bold, italics, headings, bullet points etc.\n", "- Raw NBConvert cells are rare, not used in this tutorial, but contain data that is neither rendered or executed except when exporting to PDF or the like\n", "- Heading, a soon to be deprecated feature, which gives you a heading in the Markdown format (ie a heading in a markdown cell is equivalent to a heading cell)\n", "\n", "\n", "\n", "# Outline\n", "0. [Importing Modules](#0)\n", "1. [Basic Analysis](#1)\n", " - [1.1 Basic Analysis](#11)\n", " - [1.1.1 Load Data](#111)\n", " - [1.1.2 Define some initial \"guess\" models](#112)\n", " - [1.1.3 Optimize Model](#113)\n", " - [1.2 Viterbi](#12)\n", " - [1.3 Likelihood Estimators](#13)\n", " - [1.4 Analysis Workflow](#14)\n", " - [1.5 `H2MM_arr` Function](#15)\n", "2. [The `h2mm_model` object](#2)\n", " - [2.1 Accessing core `h2mm_model` parameters](#21)\n", " - [2.2 Review of `h2mm_model` Fields](#22)\n", " - [2.3 Changing Values in a Model](#23)\n", " - [2.4 `h2mm_model` Methods](#24)\n", " - [2.4.1 The `.copy()` method](#241)\n", " - [2.4.2 The `.evaluate` Method](#242)\n", " - [2.4.3 The `.optimize` Method](#243)\n", " - [2.4.4 `EM_H2MM_C` vs `.optimize` and `H2MM_arr` vs `.evaluate`: They're basically the same](#244)\n", "3. [Format of Color and Time (Data)](#3)\n", "4. [Simulating Data](#4)\n", " - [4.1 Random Path Functions](#41)\n", " - [4.1.1 Generating a Random State Path](#411)\n", " - [4.1.2 Generating a Random Photon Path](#412)\n", " - [4.1.3 Generating States and Streams Together](#413)\n", " - [4.2 Setting the Random Seed](#42)\n", " - [4.3 Recreating Data](#43)\n", "5. [Multiparameter H2MM](#5)\n", " - [5.1 Working with more than 2 streams](#51)\n", " - [5.2 Setting up the **emission probability matrix**](#52)\n", " - [5.3 Checking a Model: `.ndet` and `.nstate` fields in `h2mm_model` objects](#53)\n", "6. [Controlling `EM_H2MM_C`](#6)\n", " - [6.1 Convergence](#61)\n", " - [6.1.1 Using `.optimize` instead](#611)\n", " - [6.2 Limits](#62)\n", " - [6.2.1 Our First Bounded Fit](#621)\n", " - [6.2.2 The `h2mm_limits` object](#622)\n", " - [6.2.3 `bounds_func` Options](#623)\n", " - [6.3 Auto-generate Initial Model with `factory_h2mm_model` function](#63)\n", " - [6.3.1 Use `factory_h2mm_model` to make initial model](#631)\n", " - [6.3.2 Bounds and `factory_h2mm_model`](#632)\n", " - [6.3.3 Distribution kwargs in `factory_h2mm_model`](#633)\n", " - [6.3.4 Full analysis with bounds and `factory_h2mm_model`](#634)\n", " - [6.4 \"Continuous\" Analysis](#64)\n", " - [6.5 While loop analysis](#65)\n", "7. [Customizing `EM_H2MM_C`](#7)\n", " - [7.1 Changing parallelization: the `num_cores` kwarg](#71)\n", " - [7.2 Printing `EM_H2MM_C` iterations](#72)\n", " - [7.2.1 Basic `print_func` options](#721)\n", " - [7.2.2 Customized display using your own function](#722)\n", " - [7.3 User Defined Limits Function](#73)\n", "8. [Apendices](#8)\n", " - [8.1 `h2mm_model` fields](#81)\n", " - [8.1.1 Statistical Fields](#811)\n", " - [8.1.2 `.niter`](#812)\n", " - [8.1.3 Convergence Info](#813)\n", " -[8.2 `EM_H2MM_C` and `.optimize` Arguments and Fields](#82)\n", " - [8.2.1 The `reset_niter` argument](#821)\n", "\n", "### Text color explanation\n", "Cells that are examples of \"complete\" analysis are preceded by headings in blue. The fast-learning/experience Python programmer might be able to piece everything togoether with just these cells. These cells will also be useful as \"skeletons\" for building analysis of other data. (Please use as a skeleton, not just copy paste)\n", "\n", "### List of \"skeleton code\" cells\n", "- [Optimizing a model](#optimization)\n", "- [Basic Workflow](#wholeanalysis1)\n", "- [Error Evaluation](#erroreval) Note: this does not include the prior optimization steps\n", "- [Analysis using auto-generated models](#fhmwrkflow) Note: this and the next blue cell are very similar\n", "- [Analysis with bounds and autogenerated models](#634)\n", "- [Continuous Analysis Worklfow](#optwrkflow)\n", "- [While loop based Workflow](#finaldemo)\n", "\n", "Other colors indicate other things\n", "- Gray indicates a review of information mostly already presented- good to review, but skippable if you've been paying attention.\n", "- Green indicates that the next cell is introducing a property, and is in some way new. \n", "- Brown indicates advice and/or warnings for researchers, especially those who are coming from Matlab.\n", "- Teal indicates this text is an explanation for experienced Python programmers, talking about some deep level behavior of the `H2MM_C` API\n", "- Magenta (of which there is only one) is a deeper explanation of the **C** level inner workings of the `H2MM_C` API\n", "\n", "This tutorial is long, bellow are the most impoartant chapters/sections. You might be able to skip the others.\n", "\n", "[Chapter 1](#1) is the most essential, and should be enough to get you off the ground- don't skip it!\n", "\n", "[Chapter 5](#5) is **very** helpful in understanding multi-parameter fitting, and how the data is changed, but does not introduce anything that couldn't be reasoned out from chapter 1.\n", "\n", "The remaining chapters are for \"advanced\" fitting. So they can wait for later, but you should look into them if you want to be able to use the `H2MM_C` API to its full potential.\n", "\n", "[Section 6.2](#62) is important reading for those who want to use mpH2MM with $\\mu$sALEX experiments.\n", "\n", "[Section 7.2](#72) and [section 7.3](#73) detail how to write custom printing, and very powerful, but also very dangerous, limits functions respectively.\n", "\n", "\n", "# 0. Importing Modules\n", "\n", "First, we need to import some modules, we're going to need:\n", "- `numpy` for some array inputs, we use `import numpy as np` so that whenever we use numpy we can call it with the prefix `np.`. This is standard practice in scientific computing.\n", "- `matplotlib` for plotting graphs so we can visualize our results, we use `from matplotlib import pyplot as plt` because it turns out that we'll only need the `pyplot` subset of `matplotlib`, and `plt` is much less wordy than typing `matplotlib.pyplot` over and over. Like how we imported `numpy`, this is standard practice in scientific computing. \n", "- `H2MM_C` the API for mpH2MM, we use `from H2MM_C import *` so all functions/classes are available specified without an addition `H2MM_C.` identifier\n", "\n", "We're also importing\n", "- `os` so we can load our sample datafile, its name is already so short, we just call it as `import os`\n", " - in real analysis, this will likely be replaced by a module like FRETbursts, which will analyze the raw data.\n", " - Often in real analysis you will still need many of the functions in `os`, so you'll probably import `os` anyways, even if you are using FRETBursts" ] }, { "cell_type": "code", "execution_count": null, "id": "civic-frequency", "metadata": {}, "outputs": [], "source": [ "#import the needed modules\n", "import os\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "from H2MM_C import *" ] }, { "cell_type": "markdown", "id": "c34310d7", "metadata": {}, "source": [ "\n", "# 1. Basic Analysis\n", "\n", "mpH2MM is fundamentally an optimization method. It's main goal is to find the \"optimal\" **model** that describes a set of **data**. The **model** assumes the system is transitioning between a finite number of **states**. It is further assumed that the **states** are \"hidden\" and that the data instead comes in the form of a set of data points, which we will call **photons** because that is what they are when we are working with single molecule confocal data, the system for which H2MM was designed. Each **photon** has a **time** and **detector** (or **stream**) associated with it. Each **state** has a certain probability of producing a photon, in each **stream**. Usually you will have many short sequences that all represent the same system, which we will refer to as **bursts**, again, because these usually represent bursts of photons caused by single molecules diffusing in and out of the confocal volume.\n", "\n", "Each **model** is composed of a set of arrays:\n", "1. The **initial probability matrix**, this is the least important array, which describes the probaility of the data beginning in each state.\n", "2. The **transition probability matrix**, the likelihood of one state transitioning to another\n", "3. The **emission probability matrix**, the likelihood of a state producing a data point in a given stream\n", "\n", "\"Optimization\" is specifically the optimization of the loglikelihood, basically, it is the likelihood that the model describes the data, so optimization is a process to find the model that has the highest, i.e. the best chance of describing the data.\n", "\n", "\n", "## 1.1 Run Analysis\n", "So, the basic process is:\n", "1. Load and/or process data\n", "2. Define an initial **model**\n", "3. Optimize **model** against **data**\n", "\n", "Since mpH2MM is usually used for confocal burst data, the **streams** are the detectors, and data is usually composed of a set of bursts of photons. So the first step is usually taking your burst selection, and converting it into a form that `H2MM_C` can accept.\n", "\n", "\n", "### 1.1.1 Load Data\n", "#### First, we need some data, let's load it from the provided file:\n", "- This is just a sample data for us to work with, you'd usually use a module like FRETBursts to process the raw data and slice it into bursts.\n", "- If you are using Matlab for your initial data analysis, then you will want to export the data to an ascii text file, and load it into Python using code similar to the cell bellow (make sure you get the format correct, especially note that python indexes from 0, while Matlab indexes from 1, so you want to make sure that photon streams are indexed as 0,1,2... and not as 1,2,3...\n", "\n", "#### The **data** is imported into two lists: `color` and `times`\n", "- `color` is a list* (it can also be a tuple) of NUMPY arrays, each array is a **burst** which contains the **stream** indices (0,1,2...) of the photons *(Note that this is not exclusively based on the **detector** as two streams may also be distinguished by the excitation period, so one **detector** might give rise to two **streams**: one for donor excitaiton, and one for acceptor excitation)*\n", "- `times` is a list* (it can also be a tuple) of NUMPY arrays, much like in `color` each burst gets an array, but instead of containing the **stream** indeces, these are the photon arrival times (arbitrary units). Make sure that the elements of `color` and `times` match, i.e. they're the same length, we'll be more specific later.\n", "\n", "*For those coming from Matlab, Python lists are very similar to Matlab Cell Arrays, Python tuples are similar, except that they cannot be changed after their creation.\n", "\n", "### Load Bursts:" ] }, { "cell_type": "code", "execution_count": null, "id": "closing-yeast", "metadata": {}, "outputs": [], "source": [ "# This cell will only be used if you export your data to an ascii file, \n", "# if you are using a software package such as FRETBursts, a different\n", "# function (for FRETBursts, look at our other Jupyter Notebooks, it is\n", "# the data_sort function)\n", "times = []\n", "color = []\n", "i = 0\n", "# load the data\n", "with open('sample_data_2det.txt') as f:\n", " for line in f: # this reads each line one at a time\n", " if i % 2 == 0: # even lines are times\n", " times.append(np.array([int(x) for x in line.split()],dtype='Q'))\n", " else: # odd lines are color\n", " color.append(np.array([int(x) for x in line.split()],dtype='L'))\n", " i += 1\n", "print(\"Done loading data\")" ] }, { "cell_type": "markdown", "id": "9df5fabc", "metadata": {}, "source": [ "\n", "### 1.1.2 Define some initial \"guess\" models\n", "\n", "In `H2MM_C` the **model** is a special object: the `h2mm_model` object (for those coming from Matlab, think of it like a structure with predefined fields, but remember it has a lot more functionality built into it than a Matlab structure).\n", "There are 3 components that define every **model,** these are specified as NUMPY array (floating point). They are all row stochastic, i.e. each row sums to 1, and negative values are not allowed:\n", "- **Initial Probability**: (`.prior`) the probability of the system being in each state at the beginning of a **burst**.\n", " - A 1D array, with length of the number of states\n", "- **Transition Probability matrix**: (`.trans`) each element `[i,j]` is the probability of the system transitioning from state `i` to state `j` in one time step (i.e. the unit of time used as 1 in the times array).\n", " - A 2D square array, with dimensions of the number of states\n", " - Usually the values on the diagonal are close to 1, and off-diagonal close to 0\n", "- **Emission Probability matrix**: (`.obs`) the probability that a photon will arise from a given photon stream, i.e. the `[i,j]` element is the probability that, given the system is in state `i`, a photon will belong to photon stream `j`.\n", " - A 2D array, with dimension 0 being the number of states, and dimension 1 being the number of photon streams\n", "\n", "*A final note: when a `h2mm_model` object is created, the `H2MM_C` API automatically \"normalizes\" these input arrays, that is it makes all rows sum to 1. It is highly recommended to use already row stochastic matrices as inputs.*\n", "\n", "First, lets make the arrays for input:\n", "#### Define Model Arrays" ] }, { "cell_type": "code", "execution_count": null, "id": "architectural-pontiac", "metadata": {}, "outputs": [], "source": [ "# first we're going to define the arrays\n", "prior = np.array([0.5,0.5]) # note: it is a 1D array, and has 2 elements, that means this model has 2 states\n", "trans = np.array([[0.999999, 1e-6],[1e-6, 0.999999]]) # note that it is 2D, and square, with 2 states\n", "obs = np.array([[0.3, 0.7],[0.5, 0.5]]) # note that it is 2D, with 2 states (rows) and 2 streams (columns)" ] }, { "cell_type": "markdown", "id": "e7772d14", "metadata": {}, "source": [ "Now we use those arrays to declare an initial model, the class of a model is `h2mm_model`, and the call to initialize an object of this type we use the format `h2mm_model(prior,trans,obs)`\n", "\n", "#### Make Model Object" ] }, { "cell_type": "code", "execution_count": null, "id": "accompanied-davis", "metadata": {}, "outputs": [], "source": [ "# make the 2 state model\n", "init_model = h2mm_model(prior,trans,obs)" ] }, { "cell_type": "markdown", "id": "690da17d", "metadata": {}, "source": [ "\n", "### 1.1.3 Optimize Model\n", "\n", "The core function of the `H2MM_C` API is called `EM_H2MM_C`, which takes 3 inputs:\n", "- The initiating `h2mm_model`, (`init_model` in this example)\n", "- The list or tuple of NUMPY arrays of photon streams, (`color` in this example)\n", "- the list or tuple of NUMPY arrays of arrival times, (`times` in this example)\n", "\n", "The function outputs a new, optimized `h2mm_model` object.\n", "\n", "The basic signature is `EM_H2MM_C(model,color,time)`.\n", "\n", "### Run the optimization" ] }, { "cell_type": "code", "execution_count": null, "id": "wooden-fitness", "metadata": {}, "outputs": [], "source": [ "# run the optimization\n", "opt_model = EM_H2MM_C(init_model,color,times)" ] }, { "cell_type": "markdown", "id": "bc3bac9d", "metadata": {}, "source": [ "### Examining the result" ] }, { "cell_type": "markdown", "id": "698e34e5", "metadata": {}, "source": [ "You can print everything out in a jupyter notebook running a cell that has only the name of the output variable: \n", "\n", "*For those who like to get deep into Python, this calls the `.__repr__()` method.*" ] }, { "cell_type": "code", "execution_count": null, "id": "imposed-wellington", "metadata": {}, "outputs": [], "source": [ "opt_model" ] }, { "cell_type": "markdown", "id": "86d7e94a", "metadata": {}, "source": [ "Now we have everything we need to run an optimization, so let's see how this looks like in one cell.\n", "\n", "\n", "### Here we see the whole analysis in one cell:" ] }, { "cell_type": "code", "execution_count": null, "id": "a6ac2bd2", "metadata": {}, "outputs": [], "source": [ "# load data\n", "times = []\n", "color = []\n", "i = 0\n", "# open the file\n", "with open('sample_data_2det.txt') as f:\n", " for line in f: # this reads each line one at a time\n", " if i % 2 == 0: # the datafiel is formate\n", " times.append(np.array([int(x) for x in line.split()],dtype='Q'))\n", " else:\n", " color.append(np.array([int(x) for x in line.split()],dtype='L'))\n", " i += 1\n", "\n", "# first we're going to define the arrays\n", "prior = np.array([0.5,0.5])\n", "trans = np.array([[0.999999, 1e-6],[1e-6,0.999999]])\n", "obs = np.array([[0.3, 0.7],[0.5, 0.5]])\n", "init_model = h2mm_model(prior,trans,obs)\n", "\n", "# optimize the model\n", "opt_model = EM_H2MM_C(init_model,color,times)\n", "# print optimized model\n", "opt_model" ] }, { "cell_type": "markdown", "id": "f4e1c06c", "metadata": {}, "source": [ "\n", "## 1.2 Viterbi\n", "\n", "Now that a model has been found, it can be useful to find the *most likely* state of each **photon** in a **burst**. This is achieved using the *Viterbi* algorithm.\n", "\n", "There are 2 levels of *Viterbi* algorithm analysis: `viterbi_path` and `viterbi_sort`.\n", "`viterbi_sort` is a wrapper around `viterbi_path` which returns more outputs, processing the data based on dwells.\n", "Both take the same core inputs as `EM_H2MM_C`, that is `viterbi_path(model,color,times)` and `viterbi_path(model,color,times)`.\n", "The basic outputs of `viterbi_path` are (in order)\n", "- *Path*: the most likely state path through the data- a list of 1D int numpy arrays\n", "- *Scale*: the posterior probability of each photon, ie, the *likelihood* of the given photon belonging to the state given in `path`- a list of 1D float numpy arrays\n", "- *Loglik* the logliklihood of each path, given as a 1D NUMPY array, with each element cooresponding to a burst.\n", "- *ICL*: the integrated complete likelihood ($ICL$, see [section 1.3.2](#132)) of the path through all the data, basically $BIC$ (see [section 1.3.1](#131)) of the *Viterbi* path \n", "\n", "The outputs of `viterbi_sort` are reordered slightly, with the ICL put first, and outputs placed after the outputs of `viterbi_path`see the docstring for more details (run `viterbi_sort?`).\n", "\n", "#### Run *Viterbi* algorithm" ] }, { "cell_type": "code", "execution_count": null, "id": "b4fe0dc5", "metadata": {}, "outputs": [], "source": [ "# run viterbi_path on optimized model\n", "path, scale, loglik, icl = viterbi_path(opt_model,color,times)" ] }, { "cell_type": "markdown", "id": "09fac253", "metadata": {}, "source": [ "Now, let's look at an element of `path` to see what it looks like" ] }, { "cell_type": "code", "execution_count": null, "id": "e1f8b357", "metadata": { "scrolled": true }, "outputs": [], "source": [ "path[0]" ] }, { "cell_type": "markdown", "id": "00cc5c24", "metadata": {}, "source": [ "Also note that each array in `path` and `scale` match the size of the equivalent arrays in `color` and `times`" ] }, { "cell_type": "code", "execution_count": null, "id": "3856d5c3", "metadata": {}, "outputs": [], "source": [ "for p, s, t in zip(path,scale,times):\n", " assert p.size == s.size == t.size\n", "print(\"All arrays match size of equivalent arrays in times\")" ] }, { "cell_type": "markdown", "id": "74cf3a89", "metadata": {}, "source": [ "\n", "## 1.3 Likelihood Estimators\n", "\n", "It was mentioned that `EM_H2MM_C` maximizes the *loglikelihood* of a model given a set of data. This is a field in the `h2mm_model` object, simply named `.loglik`.\n", "\n", "*NOTE: when an initiating model is made with `h2mm_model(prior,trans,obs)` method, the initial model has a `.loglik` of `-inf`.*" ] }, { "cell_type": "code", "execution_count": null, "id": "a941901f", "metadata": {}, "outputs": [], "source": [ "opt_model.loglik" ] }, { "cell_type": "markdown", "id": "8c72a95c", "metadata": {}, "source": [ "Remember that it was mentioned that the optimization of `EM_H2MM_C` doesn't change the number of **states**? This is important, because unless you know for sure that there are $n$ **states**, models with different numbers of **states** should be optimized and then compared.\n", "\n", "So use it is generally better to proceed as follows:\n", " 1. Define several **models** with different numbers of **states**\n", " 2. Run optimization on **data** using all **models**\n", " 3. Compare optimized **models** and select the best one\n", "\n", "How to compare the models? Don't just choose the model you \"like,\" rather, it's best to use a statistical discriminator (or several). The *loglikelihood* however, is not a good discriminator. Generally speaking, it is good only for comparing models with the same number of states. Models with more states basically always have better *loglikelihoods.* But there are ways to make working statistical discriminators. \n", "\n", "There are two that `H2MM_C` has built-in: the $BIC$ and $ICL$\n", "\n", "\n", "### 1.3.1 Bayes Information Criterion\n", "\n", "The *Bayes Information Criterion* ($BIC$), which penalizes the *loglikelihood* based on the number of free parameters (related to the number of *states*, and factoring in the number of **streams**.\n", "\n", "Here is the statistical definition of the $BIC$:\n", "\n", "$ BIC = -2 LL + K \\ln(n)$\n", "\n", "Where $LL$ is the *loglikelihood*, $K$ is the number of free parameters, and $n$ is the number of **photons** in the data set.\n", "\n", "Since the $BIC$ is based purely on the *loglikelihood* and the data against which the **model** was optimized, the $BIC$ is a field available in the `h2mm_model`: accessed through `.bic`\n", "\n", "\n", "### 1.3.2 Integrated Complete Likelihood\n", "\n", "The *Integrated Complete Likelihood* ($ICL$) is basically the same, but based on the *loglikelihood* not of the **model**, but rather of the *most likely state path* based on the *Viterbi* algorithm.\n", "\n", "$ ICL = -2 LL_{Viterbi} + K \\ln(n)$\n", "\n", "We find the $ICL$ is a better discriminator that $BIC$. Unfortunately, unlike the $BIC$, since $ICL$ is based on the *most likely state path*, it requires running `viterbi_path` or `viterbi_sort`, and is therefore not a field in the `h2mm_model` object.\n", "\n", "So, let's see how to access and calculate both of these (if you were paying close attention, you'll already know how to find the $ICL$)\n", "\n", "#### Finding the BIC:" ] }, { "cell_type": "code", "execution_count": null, "id": "525fe11e", "metadata": {}, "outputs": [], "source": [ "# the BIC:\n", "opt_model.bic" ] }, { "cell_type": "markdown", "id": "f8c5370c", "metadata": {}, "source": [ "### And the ICL\n", "(remember, we have to run *Viterbi*)" ] }, { "cell_type": "code", "execution_count": null, "id": "fa9a6937", "metadata": {}, "outputs": [], "source": [ "# calculate viterbi, the ICL is the last output\n", "path, scale, loglik, icl = viterbi_path(opt_model,color,times)\n", "# print the ICL:\n", "icl" ] }, { "cell_type": "markdown", "id": "ff6bf7b7", "metadata": {}, "source": [ "\n", "## 1.4 Analysis Workflow\n", "\n", "Now we have all the tools we need to properly analyze a set of data, and decide what **model** best describes it:\n", "1. Load and/or process data (we recommend using *FRETBursts*)\n", "2. Define several initial **models** with different numbers of states\n", "3. Optimize each **model** against the **data**\n", "4. Calculate $ICL$ with the *Viterbi* algorithm\n", "5. Compare $ICL$ (and potentially $BIC$) to select the best model, the minimum is considered best\n", "\n", "\n", "#### So, let's see this in action:" ] }, { "cell_type": "code", "execution_count": null, "id": "b8f7fdb2", "metadata": {}, "outputs": [], "source": [ "# load data (this is the same as before)\n", "times = []\n", "color = []\n", "i = 0\n", "# open the file\n", "with open('sample_data_2det.txt') as f:\n", " for line in f: # this reads each line one at a time\n", " if i % 2 == 0: # the datafiel is formate\n", " times.append(np.array([int(x) for x in line.split()],dtype='Q'))\n", " else:\n", " color.append(np.array([int(x) for x in line.split()],dtype='L'))\n", " i += 1\n", "\n", "\n", "# define several models, with different numbers of states (this is different from previous cells)\n", "# 1state model\n", "prior_1state = np.array([1.0])\n", "trans_1state = np.array([[1.0]])\n", "obs_1state = np.array([[0.5,0.5]])\n", "init_model_1state = h2mm_model(prior_1state,trans_1state,obs_1state)\n", "# 2 state model\n", "prior_2state = np.array([0.5,0.5])\n", "trans_2state = np.array([[0.999999, 1e-6],[1e-6, 0.999999]])\n", "obs_2state = np.array([[0.3, 0.7],[0.5, 0.5]])\n", "init_model_2state = h2mm_model(prior_2state,trans_2state,obs_2state)\n", "\n", "# 3 state model\n", "prior_3state = np.array([0.3, 0.3, 0.4])\n", "trans_3state = np.array([[0.999998,1e-6,1e-6],[1e-6,0.999998,1e-6],[1e-6,1e-6,0.999998]])\n", "obs_3state = np.array([[0.3,0.7],[0.5,0.5],[0.7,0.3]])\n", "init_model_3state = h2mm_model(prior_3state,trans_3state,obs_3state)\n", "\n", "# 4 state model\n", "prior_4state = np.array([0.2,0.2,0.2,0.4])\n", "trans_4state = np.array([[0.999997,1e-6,1e-6,1e-6],[1e-6,0.999997,1e-6,1e-6],[1e-6,1e-6,0.999997,1e-6],[1e-6,1e-6,1e-6,0.999997]])\n", "obs_4state = np.array([[0.2,0.8],[0.4,0.6],[0.6,0.4],[0.8,0.2]])\n", "init_model_4state = h2mm_model(prior_4state,trans_4state,obs_4state)\n", "\n", "\n", "# Optimize each model\n", "opt_model_1state = EM_H2MM_C(init_model_1state,color,times)\n", "opt_model_2state = EM_H2MM_C(init_model_2state,color,times)\n", "opt_model_3state = EM_H2MM_C(init_model_3state,color,times)\n", "opt_model_4state = EM_H2MM_C(init_model_4state,color,times)\n", "\n", "\n", "# calculate ICL\n", "path_1state, scale_1state, ll_1state, icl_1state = viterbi_path(opt_model_1state,color,times)\n", "path_2state, scale_2state, ll_2state, icl_2state = viterbi_path(opt_model_2state,color,times)\n", "path_3state, scale_3state, ll_3state, icl_3state = viterbi_path(opt_model_3state,color,times)\n", "path_4state, scale_4state, ll_4state, icl_4state = viterbi_path(opt_model_4state,color,times)" ] }, { "cell_type": "markdown", "id": "79ed82d8", "metadata": {}, "source": [ "Now, let's examine the results, first, let's look at the $BIC$:" ] }, { "cell_type": "code", "execution_count": null, "id": "cb892c7e", "metadata": {}, "outputs": [], "source": [ "print(f\"1 state BIC = {opt_model_1state.bic}\")\n", "print(f\"2 state BIC = {opt_model_2state.bic}\")\n", "print(f\"3 state BIC = {opt_model_3state.bic}\")\n", "print(f\"4 state BIC = {opt_model_4state.bic}\")" ] }, { "cell_type": "markdown", "id": "3269b3ff", "metadata": {}, "source": [ "And then the $ICL$:" ] }, { "cell_type": "code", "execution_count": null, "id": "9a886b2b", "metadata": {}, "outputs": [], "source": [ "print(f\"1 state ICL = {icl_1state}\")\n", "print(f\"2 state ICL = {icl_2state}\")\n", "print(f\"3 state ICL = {icl_3state}\")\n", "print(f\"4 state ICL = {icl_4state}\")" ] }, { "cell_type": "markdown", "id": "dfe70add", "metadata": {}, "source": [ "Let's use Matplolib to graphically represent this:" ] }, { "cell_type": "code", "execution_count": null, "id": "690cc989", "metadata": {}, "outputs": [], "source": [ "plt.scatter([1,2,3,4],[opt_model_1state.bic, opt_model_2state.bic, opt_model_3state.bic, opt_model_4state.bic])\n", "plt.title(\"BIC\")" ] }, { "cell_type": "code", "execution_count": null, "id": "4b8d6aad", "metadata": {}, "outputs": [], "source": [ "plt.scatter([1,2,3,4],[icl_1state, icl_2state, icl_3state, icl_4state])\n", "plt.title(\"ICL\")" ] }, { "cell_type": "markdown", "id": "54880274", "metadata": {}, "source": [ "The $ICL$ of the 2 state model is at a minimum, therefore it is most likely the ideal model. \n", "\n", "The $BIC$ continuously improves, but in our testing, this is usually the case, however, the improvement decreases with successive states- so you might choose a threshold, or like *Lerner2018*, use a modified $BIC$.\n", "\n", "*We've also noticed that the principle of the more data the better is especially true for H2MM analysis. So if you're seeing plateaus, especially in $ICL$, consider acquiring more data to make things clear.*\n", "\n", "So, let's see the ideal model:" ] }, { "cell_type": "code", "execution_count": null, "id": "439ee1c4", "metadata": {}, "outputs": [], "source": [ "opt_model_2state" ] }, { "cell_type": "markdown", "id": "ffcc6e1b", "metadata": {}, "source": [ "\n", "## 1.5 `H2MM_arr` Function\n", "\n", "Sometimes you don't want to optimize a model, but rather just get the *loglikelihood* or $BIC$ out of that model. Maybe you want to get a bunch of loglikelihoods against the same set of data.\n", "\n", "For this, there is the `H2MM_arr` function, which has a very similar signature to that of `EM_H2MM_C`, `viterbi_path`, and `viterbi_path`: `H2MM_arr(models,color,time)` The main difference is that `models` can be a `list`, `tuple`, NUMPY array of `h2mm_model` objects, or just a `h2mm_model` object.\n", "\n", "#### This can be used as a way to see how sharply the *loglikelihood* peaks around it's optimal value, which can help you estimate error bars.\n", "\n", "To do this, we'll first extract the model parameters, (`.prior`, `.trans`, and `.obs`) into separate arrays, and make several initial **models** from those, each with one value slightly tweaked (and renormalized), and then evaluate the *loglikelihood* with `H2MM_arr`:\n", "\n", "#### Setup the new models" ] }, { "cell_type": "code", "execution_count": null, "id": "a9d8d818", "metadata": {}, "outputs": [], "source": [ "# extract prior, trans and obs\n", "prior_opt = opt_model_2state.prior\n", "trans_opt = opt_model_2state.trans\n", "obs_opt = opt_model_2state.obs\n", "\n", "# make copies of the array so the tweak can be applied\n", "# we need to do \"+ 0\" so that the tweak arrays are separate copies in memory, you could also use the .copy() method\n", "trans_tweak1 = trans_opt + 0\n", "trans_tweak2 = trans_opt + 0\n", "trans_tweak3 = trans_opt + 0\n", "trans_tweak4 = trans_opt + 0\n", "\n", "# let's tweak the value of trans[0,1], lets adjust it by 1e-7 up and down\n", "trans_tweak1[0,1] = trans_opt[0,1] - 2e-7\n", "trans_tweak2[0,1] = trans_opt[0,1] - 1e-7\n", "trans_tweak3[0,1] = trans_opt[0,1] + 1e-7\n", "trans_tweak4[0,1] = trans_opt[0,1] + 2e-7\n", "\n", "# to renormalize, we'll have to do the oposite to [0,0]\n", "trans_tweak1[0,0] = trans_opt[0,0] + 2e-7\n", "trans_tweak2[0,0] = trans_opt[0,0] + 1e-7\n", "trans_tweak3[0,0] = trans_opt[0,0] - 1e-7\n", "trans_tweak4[0,0] = trans_opt[0,0] - 2e-7\n", "\n", "# make the models\n", "model_tweak1 = h2mm_model(prior_opt,trans_tweak1,obs_opt)\n", "model_tweak2 = h2mm_model(prior_opt,trans_tweak2,obs_opt)\n", "model_tweak3 = h2mm_model(prior_opt,trans_tweak3,obs_opt)\n", "model_tweak4 = h2mm_model(prior_opt,trans_tweak4,obs_opt)" ] }, { "cell_type": "markdown", "id": "abd9ac5a", "metadata": {}, "source": [ "#### Now run `H2MM_arr` (Error Evaluation)" ] }, { "cell_type": "code", "execution_count": null, "id": "7954785c", "metadata": {}, "outputs": [], "source": [ "# put the new models into a list\n", "model_list = [model_tweak1, model_tweak2, opt_model_2state, model_tweak3, model_tweak4]\n", "out_list = H2MM_arr(model_list,color,times)" ] }, { "cell_type": "markdown", "id": "110cec29", "metadata": {}, "source": [ "#### Examine the results:\n", "*NOTE: the cell bellow takes advantage of feature of Python called list comprehension, if the syntax doens't make sense, check out this tutorial. it's a really powerful feature*." ] }, { "cell_type": "code", "execution_count": null, "id": "5f54cc49", "metadata": {}, "outputs": [], "source": [ "loglik = [model.loglik for model in out_list]\n", "trans_rate = [model.trans[0,1] for model in out_list]\n", "plt.scatter(trans_rate,loglik)" ] }, { "cell_type": "markdown", "id": "e1107738", "metadata": {}, "source": [ "\n", "# 2. The `h2mm_model` object\n", "\n", "\n", "## 2.1 Accessing core `h2mm_model` parameters\n", "\n", "We've already seen how to look at the *loglikelihood* and $BIC$ from an optimized `h2mm_model` with `model.loglik` and `model.bic`. But you can also look at the arrays that make up the model- those all important **initial probability**, **transition probability** and **emission probability** matrices. They are named basically the same as how we defined them:\n", " - `.prior` for the **initial probability**\n", " - `.trans` for the **transition probability**\n", " - `.obs` for the **emission probability**\n", "\n", "Whenever you call one of these properties, Python will return a NUMPY array. You can even assign it to a new variable, but note that if you do, changing the values in that variable will **not** change the values in the model, this is because under the hood, Python is actually handing you a copy of the array.\n", "\n", "#### The initial model" ] }, { "cell_type": "code", "execution_count": null, "id": "lightweight-nickel", "metadata": {}, "outputs": [], "source": [ "prior_in = init_model_2state.prior\n", "prior_in" ] }, { "cell_type": "code", "execution_count": null, "id": "southwest-coral", "metadata": {}, "outputs": [], "source": [ "trans_in = init_model_2state.trans\n", "trans_in" ] }, { "cell_type": "code", "execution_count": null, "id": "fleet-springer", "metadata": {}, "outputs": [], "source": [ "obs_in = init_model_2state.obs\n", "obs_in" ] }, { "cell_type": "markdown", "id": "lucky-heritage", "metadata": {}, "source": [ "#### The optimized model" ] }, { "cell_type": "code", "execution_count": null, "id": "interim-direction", "metadata": {}, "outputs": [], "source": [ "prior_opt = opt_model_2state.prior\n", "prior_opt" ] }, { "cell_type": "code", "execution_count": null, "id": "accompanied-tonight", "metadata": {}, "outputs": [], "source": [ "trans_opt = opt_model_2state.trans\n", "trans_opt" ] }, { "cell_type": "code", "execution_count": null, "id": "focal-independence", "metadata": {}, "outputs": [], "source": [ "obs_opt = opt_model_2state.obs\n", "obs_opt" ] }, { "cell_type": "markdown", "id": "789f3c15", "metadata": {}, "source": [ "Now, we can change values in these arrays, or do whatever else you want. Since these are **copies** we will not be changing anything in the **models**.\n", "\n", "One of the most obvious things to do is to convert the the **transition probability matrix** into units of $s^{-1}$, instead of whatever the clock of the measurement is.\n", "*Remember that `H2MM_C` treats the arrival times as integers, so you have to know what the difference in time is between time 1 and time 2. In many nsALEX experiments, the file will store arrival times this way, and will include a clock rate, the time between 1 and 2 is so many nanoseconds.*\n", "\n", "Let's choose a value common to nsALEX measurements: $50ns$, in seconds this is $5.0*10^{-5}s$" ] }, { "cell_type": "code", "execution_count": null, "id": "cdae1d51", "metadata": {}, "outputs": [], "source": [ "# suppose the resoluiton of times is 50e-6 (50MHz), then dividing by that will give the transition rate in s^-1\n", "clk_p = 5.0e-5\n", "trans_opt = trans_opt/clk_p\n", "\n", "print(\"trans_out = \")\n", "print(trans_opt)\n", "print(\"h_out.trans = \")\n", "print(opt_model_2state.trans)" ] }, { "cell_type": "code", "execution_count": null, "id": "6e585409", "metadata": { "scrolled": true }, "outputs": [], "source": [ "if np.all(trans_opt == opt_model_2state.trans):\n", " print(\"Changing trans_out changed opt_model_2state.trans\")\n", "else:\n", " print(\"Changing trans_out left opt_model_2state.trans the same\")" ] }, { "cell_type": "markdown", "id": "9d86132f", "metadata": {}, "source": [ "\n", "## 2.2 Review of `h2mm_model` Fields\n", "\n", "Now let's review what we've learned about the `h2mm_model` fields:\n", "1. Model Definition: these are the values that define a model:\n", " - `.prior` the **initial probability matrix** is a 1D array, and has length of the number of states in the model, the sum of all values in this matrix should be 1, and none negative- it represents the probability of the system being in a given state at the beginning of a burst, so the system will be in a state at the beginning, and negative probabilities are meaningless)\n", " - `.trans` the **transition probability matrix** a 2D array, both dimensions being the number of states. Each row should sum to 1, as it represents the probability that if the system is in a state, what is the probability to transition to another state, so the system must transition to a state, even if it is to the same state (in fact staying in the same state is the most likely, which is why the diagonal of this matrix is usually close to 1)\n", " - `.obs` the **emission probability matrix** a 2D array with rows as the number of **states**, and columns the number of **streams** (see chapter 3, so far we have only dealt with 2 streams), the rows should sum to one, and none negative, as each row is the probability that if the system is in that state, it will emit a photon.\n", "2. Likelihood Estimators (calculated from **data**, `.prior`, `.trans`, and `.obs`):\n", " - `.loglik`the *loglikelihood* (log of likelihood), of the model, useful for comparing two models with the same number of **states**, run against the same set of data. *Initial models (before calculation of loglik against **data**) have a value of `-inf`*\n", " - `.bic` the Bayesian Information Criterion, an assessment of the likelihood of the model, with penalties for the number of **states**, allows comparison of models with different numbers of **states**, (run against the same set of data). *Initial models (before calculation of loglik against **data** will raise an error if this field is called.*\n", "\n", "\n", "## 2.3 Changing Values in a Model\n", "\n", "When you access the values in an `h2mm_model` object, you do not change them.\n", "\n", "But what if you have a new set of data, or something else, and you would like to change those values? Well, that can be done by passing an array into the field: \n", "\n", "- `model.prior = new_prior`\n", "- `model.trans = new_trans`\n", "- `model.obs = new_obs`\n", "\n", "*NOTE the new array must be the **same** size, otherwise an error will be raised.*\n", "\n", "It should also be noted that the `H2MM_C` API automatically resets the `.loglik` to `-inf` and the `.bic` will raise a warning, when you set the value of one of these fields. This is good because the model is no longer the same, and therefore a new *loglikelihood* must be calculated. Essentially when one of these fields is set, the model should be treated as an initial model, and this is how the `H2MM_C` API treats it.\n", "\n", "#### Change the `.trans` property of an `h2mm_model` object:" ] }, { "cell_type": "code", "execution_count": null, "id": "f50e9aac", "metadata": {}, "outputs": [], "source": [ "# make a new trans array\n", "trans_new = np.array([[0.9999999,1e-7],[2e-7,0.9999998]])\n", "# assign to the model\n", "opt_model_2state.trans = trans_new\n", "# display the model\n", "opt_model_2state" ] }, { "cell_type": "markdown", "id": "2f3ea5dc", "metadata": {}, "source": [ "Note that the `.loglik` is now reset to `-inf`:" ] }, { "cell_type": "code", "execution_count": null, "id": "d4573a56", "metadata": {}, "outputs": [], "source": [ "opt_model_2state.loglik" ] }, { "cell_type": "markdown", "id": "3972b0fe", "metadata": {}, "source": [ "So, we should re-optimize this model.\n", "\n", "\n", "Since the model is no longer optimized, we should run `H2MM_arr` to get the new *loglikelihood.*" ] }, { "cell_type": "code", "execution_count": null, "id": "71806dcf", "metadata": {}, "outputs": [], "source": [ "opt_model_2state = H2MM_arr(opt_model_2state,color,times)\n", "opt_model_2state" ] }, { "cell_type": "markdown", "id": "85711587", "metadata": {}, "source": [ "Now the *loglikelihood* is real:" ] }, { "cell_type": "code", "execution_count": null, "id": "25aab7b6", "metadata": {}, "outputs": [], "source": [ "opt_model_2state.loglik" ] }, { "cell_type": "markdown", "id": "c2fdf0a0", "metadata": {}, "source": [ "\n", "## 2.4 `h2mm_model` Methods\n", "\n", "To Matlab users: one of the biggest differences between Python and Matlab is that Python is \"object oriented.\" Practically this means you can execute functions on a variable with the dot (\"`.`\") notation, which the functions available depend on the type of the variable, and these are generally called \"methods\"\n", "\n", "The `h2mm_model` object has a few methods, let's have a look at how to use them.\n", "\n", "\n", "### 2.4.1 The `.copy()` method\n", "\n", "In python, if we just execute `new_model = old_model` where `old_model` is a `h2mm_model` object, all we have really done is let ourselves refer to `old_model` in a new way. This means that if we change the `.prior`, `.trans` or `.obs` arrays in `new_model`, we will change those values in `old_model`.\n", "\n", "#### Demonstrate Problem of using `new_model = old_model`" ] }, { "cell_type": "code", "execution_count": null, "id": "ee49d2d1", "metadata": {}, "outputs": [], "source": [ "# get the current state of the trans model\n", "old_trans = opt_model_2state.trans\n", "# use the standard = operator without copy()\n", "new_model = opt_model_2state\n", "new_trans = np.array([[0.9999995,5e-7],[5e-7,0.9999995]])\n", "new_model.trans = new_trans\n", "\n", "# display the old model trans matrix, note that it changed\n", "print(\"The current state of opt_model_2state.trans = \")\n", "print(opt_model_2state.trans)\n", "print(\"The old state of opt_model_2state.trans = \")\n", "print(old_trans)\n", "if np.all(old_trans == opt_model_2state.trans):\n", " print(\"The trans array is unchanged\")\n", "else:\n", " print(\"The trans array is changed\")" ] }, { "cell_type": "markdown", "id": "a2b30b15", "metadata": {}, "source": [ "Python has the `is` operator, which tells us if two variables point to the same object. Let's use it to prove the point a litle more clearly." ] }, { "cell_type": "code", "execution_count": null, "id": "d457df2e", "metadata": {}, "outputs": [], "source": [ "new_model is opt_model_2state" ] }, { "cell_type": "markdown", "id": "7b771c70", "metadata": {}, "source": [ "The `H2MM_C` API includes the `.copy()` method, which returns a **copy** of the entire model, with its values stored separately in memory (for the experienced Python programmer, it makes a *deep copy* of the data)." ] }, { "cell_type": "code", "execution_count": null, "id": "6815cc61", "metadata": {}, "outputs": [], "source": [ "opt_model_2state = EM_H2MM_C(opt_model_2state,color,times)" ] }, { "cell_type": "markdown", "id": "0ae0a1b1", "metadata": {}, "source": [ "#### Use the `.copy()` method" ] }, { "cell_type": "code", "execution_count": null, "id": "5aea2bb2", "metadata": {}, "outputs": [], "source": [ "# get the current state of the trans model\n", "old_trans = opt_model_2state.trans\n", "# use the copy() method\n", "new_model = opt_model_2state.copy()\n", "new_trans = np.array([[0.9999995,5e-7],[5e-7,0.9999995]])\n", "new_model.trans = new_trans\n", "\n", "# display the old model trans matrix, note that it changed\n", "print(\"The current state of opt_model_2state.trans = \")\n", "print(opt_model_2state.trans)\n", "print(\"The old state of opt_model_2state.trans = \")\n", "print(old_trans)\n", "if np.all(old_trans == opt_model_2state.trans):\n", " print(\"The trans array is unchanged\")\n", "else:\n", " print(\"The trans array is changed\")" ] }, { "cell_type": "markdown", "id": "4cda9480", "metadata": {}, "source": [ "Let's use the `is` operator again:" ] }, { "cell_type": "code", "execution_count": null, "id": "678c52f7", "metadata": {}, "outputs": [], "source": [ "new_model is opt_model_2state" ] }, { "cell_type": "markdown", "id": "26b7853e", "metadata": {}, "source": [ "#### Practical example\n", "\n", "Remember when we were evaluating the error, how we had to make a bunch of arrays and make new `h2mm_model` objects using them? Well, the `.copy()` method allows us to do this without having to extract all the arrays.\n", "\n", "\n", "#### Error Evaluation with `.copy()`" ] }, { "cell_type": "code", "execution_count": null, "id": "9f874c44", "metadata": {}, "outputs": [], "source": [ "# make copies of the model for each tweaked model\n", "model_tweak1 = opt_model_2state.copy()\n", "model_tweak2 = opt_model_2state.copy()\n", "model_tweak3 = opt_model_2state.copy()\n", "model_tweak4 = opt_model_2state.copy()\n", "\n", "# makes arrays for each tweak to adjust\n", "trans_tweak1 = opt_model_2state.trans\n", "trans_tweak2 = opt_model_2state.trans\n", "trans_tweak3 = opt_model_2state.trans\n", "trans_tweak4 = opt_model_2state.trans\n", "\n", "# adjust the [0,1] element of the arrays\n", "trans_tweak1[0,1] = trans_tweak1[0,1] - 2e-7\n", "trans_tweak2[0,1] = trans_tweak2[0,1] - 1e-7\n", "trans_tweak3[0,1] = trans_tweak3[0,1] + 1e-7\n", "trans_tweak4[0,1] = trans_tweak4[0,1] + 2e-7\n", "\n", "# adjust the [0,0] so that the arrays are row stochastic\n", "trans_tweak1[0,0] = trans_tweak1[0,0] + 2e-7\n", "trans_tweak2[0,0] = trans_tweak2[0,0] + 1e-7\n", "trans_tweak3[0,0] = trans_tweak3[0,0] - 1e-7\n", "trans_tweak4[0,0] = trans_tweak4[0,0] - 2e-7\n", "\n", "# change the values in the trans arrays in the copies\n", "model_tweak1.trans = trans_tweak1\n", "model_tweak2.trans = trans_tweak2\n", "model_tweak3.trans = trans_tweak3\n", "model_tweak4.trans = trans_tweak4\n", "\n", "# make a list out of the models\n", "model_list = [model_tweak1, model_tweak2, opt_model_2state, model_tweak3, model_tweak4]\n", "# evaluate the loglik for all the tweaked models\n", "out_list = H2MM_arr(model_list,color,times)\n", "\n", "loglik = [model.loglik for model in out_list]\n", "trans_rate = [model.trans[0,1] for model in out_list]\n", "\n", "plt.scatter(trans_rate,loglik)" ] }, { "cell_type": "code", "execution_count": null, "id": "b31acf77", "metadata": {}, "outputs": [], "source": [ "trans_rate, loglik" ] }, { "cell_type": "markdown", "id": "23c5722d", "metadata": {}, "source": [ "\n", "### 2.4.2 The `.evaluate` Method\n", "\n", "Sometimes we'd like to just evaluate the *loglikelihood* of a given `h2mm_model` object against a set of **data**, and not have to make a separate object. We could use `model = H2MM_arr(model,color,time)`, but there is a cleaner way: the `.evaluate` method, which under the hood calls `H2MM_arr`, but `.evaluate` is cleaner. (This is a better way to do what we did [here](#setandeval).)\n", "\n", "*For Matlab users: this will be a good example of how methods work in object oriented programming*\n", "\n", "The signature is basically the same: `model.evaluate(color,time)`\n", "\n", "To demonstrate, let's make another model in our curve seeing how sharp the *loglikelihood* peak is around the transition rate.\n", "\n", "#### Make new tweaked model" ] }, { "cell_type": "code", "execution_count": null, "id": "1d68e965", "metadata": {}, "outputs": [], "source": [ "# make the new trans array\n", "trans_tweak5 = opt_model_2state.trans\n", "\n", "trans_tweak5[0,1] = trans_tweak5[0,1] + 3e-7\n", "trans_tweak5[0,0] = trans_tweak5[0,0] - 3e-7\n", "# update the model\n", "model_tweak5 = opt_model_2state.copy()\n", "model_tweak5.trans = trans_tweak5" ] }, { "cell_type": "markdown", "id": "5391a3a8", "metadata": {}, "source": [ "#### Get the *loglikelihood* of a model with `.evaluate`" ] }, { "cell_type": "code", "execution_count": null, "id": "335ca2c5", "metadata": {}, "outputs": [], "source": [ "model_tweak5.evaluate(color,times)\n", "model_tweak4" ] }, { "cell_type": "markdown", "id": "397d3e84", "metadata": {}, "source": [ "Look at the `.loglik` specifically:" ] }, { "cell_type": "code", "execution_count": null, "id": "55eeaf51", "metadata": {}, "outputs": [], "source": [ "model_tweak5.loglik" ] }, { "cell_type": "markdown", "id": "9d5c3c6c", "metadata": {}, "source": [ "#### Let's add the new model to the list of models\n", "Matlab users: this is a good demonstration how Python list syntax works." ] }, { "cell_type": "code", "execution_count": null, "id": "220d7333", "metadata": {}, "outputs": [], "source": [ "out_list = out_list + [model_tweak5]" ] }, { "cell_type": "markdown", "id": "5939235f", "metadata": {}, "source": [ "#### Now let's plot the updated list" ] }, { "cell_type": "code", "execution_count": null, "id": "5c69e5f9", "metadata": {}, "outputs": [], "source": [ "loglik = [model.loglik for model in out_list]\n", "trans_rate = [model.trans[0,1] for model in out_list]\n", "plt.scatter(trans_rate,loglik)" ] }, { "cell_type": "markdown", "id": "95ee12c8", "metadata": {}, "source": [ "\n", "### 2.4.3 The `.optimize` Method\n", "\n", "The `.optimize` does the same thing as `EM_H2MM_C`: find the optimal model for some **data**, but just as a method of an `h2mm_model` object. The function call is `model.optimize(color,time)` (very similar to `.evaluate`). \n", "\n", "Calling `model.optimize(color,time)` is the same as calling `model = EM_H2MM_C(model,color,time)`.\n", "\n", "So let's try this, we'll make an un-optimized model, and then optimize it against our data using `.optimize` instead of `EM_H2MM_C`.\n", "\n", "#### Make an initial model" ] }, { "cell_type": "code", "execution_count": null, "id": "80409070", "metadata": {}, "outputs": [], "source": [ "prior = np.array([0.5,0.5])\n", "trans = np.array([[0.999994,6e-6],[4e-6,0.999996]])\n", "obs = np.array([[0.8,0.2],[0.5,0.5]])\n", "model_2state = h2mm_model(prior,trans,obs)" ] }, { "cell_type": "markdown", "id": "a6315b52", "metadata": {}, "source": [ "#### Optimize with `.optimize`" ] }, { "cell_type": "code", "execution_count": null, "id": "2793b17e", "metadata": { "scrolled": true }, "outputs": [], "source": [ "model_2state.optimize(color,times)\n", "model_2state" ] }, { "cell_type": "markdown", "id": "4c473bc4", "metadata": {}, "source": [ "\n", "### 2.4.4 `EM_H2MM_C` vs `.optimize` and `H2MM_arr` vs `.evaluate`: They're basically the same\n", "\n", "The `h2mm_model` object method `.optimize` has basically the same function as the function `EM_H2MM_C`. In fact, the API uses `EM_H2MM_C` to evaluate `.optimize`. The same is true of the `.evaluate` and `H2MM_arr` function. So what you can do with one, you can do with the other. Which one you use is mostly a matter of preference. \n", "\n", "The main difference is that when you use `EM_H2MM_C`, a new `h2mm_model` object is created, so it uses a **tiny** bit more memory, but you keep the old model around in case you want to use it again. On the other hand when you use `.optimize`, the results are returned to the original model, so you save a **tiny** amount of memory, but that original model is no longer available.\n", "\n", "Since `.evaluate` is also a wrapper around `.H2MM_arr`, the same principles apply, however, since the model parameters are not changed, it might seem like using `.evaluate` is the better choice unless we are wanting to compare the same model against differet **data** sets.\n", "\n", "However, one more thing to be aware of is that the `H2MM_arr` function can also accept lists, tuples, and even NUMPY arrays of `h2mm_model` objects, and run them all against the same set of data. So it's a little more flexible in it's usage. Further, there is some preprocessing of the data that takes place before the main optimization. So if you have a lot of `h2mm_model` objects to evaluate, `H2MM_arr` will probably slightly faster, as this preprocessing takes place once if you use `H2MM_arr`, while if you `.evaluate` each model, the preprocessing will be unnecessarily repeated each for each model." ] }, { "cell_type": "markdown", "id": "needed-daughter", "metadata": {}, "source": [ "\n", "# 3. Format of Color and Time (Data)\n", "\n", "**This section is rather laborious and pedantic, many will be able to skip this section without much trouble**\n", "\n", "OK, now that we have a basic idea of how the analysis works, let's get into the details of what everything is.\n", "\n", "First, lets examine the **data**, i.e. the `color` and `times` lists (they can also be tuples).\n", "\n", "The cell bellow that runs the Python function `type` on both to get the type of both." ] }, { "cell_type": "code", "execution_count": null, "id": "occupational-laugh", "metadata": {}, "outputs": [], "source": [ "type(color), type(times)" ] }, { "cell_type": "markdown", "id": "advised-talent", "metadata": {}, "source": [ "But we should check things a little further---\n", "\n", "Think of these lists as matched pairs.\n", "First, they have to be the same length, which we can check by calling the Python function `len`:" ] }, { "cell_type": "code", "execution_count": null, "id": "endless-violence", "metadata": {}, "outputs": [], "source": [ "len(color), len(times)" ] }, { "cell_type": "markdown", "id": "secure-profile", "metadata": {}, "source": [ "Now, let's look at just one burst:" ] }, { "cell_type": "code", "execution_count": null, "id": "reported-coordinate", "metadata": {}, "outputs": [], "source": [ "color[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "inside-fossil", "metadata": {}, "outputs": [], "source": [ "times[0]" ] }, { "cell_type": "markdown", "id": "received-brazil", "metadata": {}, "source": [ "See how the `colors[0]` is switching between 0 and 1, while all the elements of `times[0]`\n", "\n", "Some things to note about bursts:\n", " - `colors[0]` and `times[0]` are NUMPY **integer** arrays, and all values are **positive** *these are even converted into unsigned integers in the code*\n", " - Values in `times[0]` are always (monotonically) increasing, this is because the represent the arrival times of the photons, so photons must be \"in order.\"\n", " - `colors[0]` is a seemingly random set of 0's and 1's. This is because they represent the index of the **stream** of each photon (i.e. 0 for Donor photons, and 1 for Acceptor photons)\n", " - The number of elements in `colors[0]` and `times[0]` is the same, because these are different aspects of the same photon\n", "\n", "Now let's demonstrate that the same is true for all elements of the `color` and `times` lists, we'll use a nifty statement called `assert`, if it evaluates as false, then Python will raise an error and stop execution:\n", "\n", "*For those coming from Matlab, the square brackets `[]` in Python indicate specifying an index, think of curly brackets and parenthesis `{}` and `()`*" ] }, { "cell_type": "code", "execution_count": null, "id": "eastern-election", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# loop over the length of the array\n", "for i in range(0,len(color)):\n", " # the type of c and t (elements of color and times) is the same, and the same as the first element of each\n", " assert type(color[i]) == type(times[i]) == type(color[0]) == type(times[0]), f\"Types of index {i} do not match\"\n", " # show that the sizes of color[i] and times[i] are the same\n", " assert color[i].size == times[i].size, f\"Sizes at index {i} do not match\"\n", "print(\"All arrays of same type and sizes of color/times match\")" ] }, { "cell_type": "markdown", "id": "80e3553b", "metadata": {}, "source": [ "Note however that different bursts do not have to be the same size:" ] }, { "cell_type": "code", "execution_count": null, "id": "55871cd9", "metadata": {}, "outputs": [], "source": [ "color[0].size, color[1].size" ] }, { "cell_type": "markdown", "id": "1a7daf3d", "metadata": {}, "source": [ "\n", "# 4. Simulating Data\n", "\n", "There is a final set of core functions in the H2MM_C API: the simulated trajectory functions.\n", "\n", "The purpose of these functions is to provide a lightweight means to generate \"simulated\" data, by which we mean generating a random set of **states** and/or **streams** based on an input **model** and set of **times**. The reason an input set of **times** is nearly always needed is that a better set of times is better generated using a too like PyBroMo. In our examples, we will instead use the random number generators functions in NUMPY to generate quick-and-dirty **time** distributions. This works only as a first approximation.\n", "\n", "\n", "## 4.1 Random Path Functions\n", "\n", "The basic process is to work backwards, taking the **model** to generate a set of **states** (like those found in *Viterbi* analysis, and then generate the **streams** from those **states**.\n", "\n", "\n", "### 4.1.1 Generating a Random State Path\n", "\n", "The function to generate a **sate path** is `sim_sparsestatepath`, and the function call is fairly simple: `sim_sparsestatepath(model,times)` where `model` is the and `h2mm_model` object, and `times` is a 1D monotonically increasing positive integer NUMPY array- the times. For this example, we will generate `times` using the NUMPY function `np.random.exponential` which generates exponentially distributed random numbers, which is to a first approximation the way photons in a burst are distributed.\n", "\n", "So lets generate a set of 50 **times**:" ] }, { "cell_type": "code", "execution_count": null, "id": "7efd2e68", "metadata": {}, "outputs": [], "source": [ "time_sim = np.cumsum(np.random.exponential(100,size=(50)))" ] }, { "cell_type": "markdown", "id": "9ca7fc73", "metadata": {}, "source": [ "#### Now we can make a random trajectory (we'll use the model from the 2 state fitting, `opt_model`:" ] }, { "cell_type": "code", "execution_count": null, "id": "2ab29975", "metadata": {}, "outputs": [], "source": [ "states = sim_sparsestatepath(opt_model,time_sim)\n", "states" ] }, { "cell_type": "markdown", "id": "67be40aa", "metadata": {}, "source": [ "Notice how the numbers are mostly consecutive? This is because, in this case, the transition rate is slow, so we expect changes of state only rarely.\n", "\n", "\n", "### 4.1.2 Generating a Random Photon Path\n", "\n", "But of course, we'd like to see what the actual data of the model might look like, as in the **streams** of the photons, not their **states**. To do this, we feed the resultant **states** into anotehr function: `sim_phtraj_from_state`, which has a very similar signature to the state generator function: `sim_phtraj_from_state(model, states)`, notice that we have swapped out the `times` array for the `states` array. This is becuase now we know the **state** of each simulated **photon**, so now we just need to use the randomly assign each photon based on the probability of a photon arrising in a given **stream** according to the **emmission probability** (`.obs`) matrix.\n", "\n", "#### Make a **stream** trajectory from a **state path**" ] }, { "cell_type": "code", "execution_count": null, "id": "185f4218", "metadata": {}, "outputs": [], "source": [ "streams = sim_phtraj_from_state(opt_model,states)\n", "streams" ] }, { "cell_type": "markdown", "id": "4179d884", "metadata": {}, "source": [ "Notice now that the numbers are much more random, this is because now we have the different **streams** all overlayed on top of each other.\n", "\n", "### 4.1.3 Generating States and Streams Together\n", "\n", "H2MM_C also provides a function that will do both steps in one, `sim_phtraj_from_times`, it's function call `sim_phtraj_from_times(model, times)`, note that this is the same as that of `sim_sparsestatepath`. But importantly, it returns both the `states` and `streams` arrays in the first and second registers.\n", "\n", "#### Generating a state path and streams from times" ] }, { "cell_type": "code", "execution_count": null, "id": "035de37a", "metadata": {}, "outputs": [], "source": [ "states, streams = sim_phtraj_from_times(opt_model,time_sim)\n", "states, streams" ] }, { "cell_type": "markdown", "id": "3394e157", "metadata": {}, "source": [ "\n", "## 4.2 Setting the Random Seed\n", "\n", "Since these simulations are based on random number generators, results will be different each time. However, if repeatability is desired, the seed of the random number generator can be set with the keyword argument `seed=int`. The same syntax is used across all three simulation functions.\n", "\n", "Bellow is an example of a use of his keyword argument:" ] }, { "cell_type": "code", "execution_count": null, "id": "d73b57fd", "metadata": {}, "outputs": [], "source": [ "states, streams = sim_phtraj_from_times(opt_model,time_sim,seed=487)\n", "states, streams" ] }, { "cell_type": "markdown", "id": "cdbbed1b", "metadata": {}, "source": [ "Now, if we use the same call, we will get the same result:" ] }, { "cell_type": "code", "execution_count": null, "id": "407b82da", "metadata": {}, "outputs": [], "source": [ "states, streams = sim_phtraj_from_times(opt_model,time_sim,seed=487)\n", "states, streams" ] }, { "cell_type": "markdown", "id": "ff41511e", "metadata": {}, "source": [ "Note: once the seed is set, it is persitent, so once it is set, so long as the inputs are exactly the same, the subsequent results will also be the same:" ] }, { "cell_type": "code", "execution_count": null, "id": "3e14108a", "metadata": {}, "outputs": [], "source": [ "# set the random seed\n", "states, streams = sim_phtraj_from_times(opt_model,time_sim,seed=487)\n", "# generate a second path\n", "states, streams = sim_phtraj_from_times(opt_model,time_sim)\n", "states, streams" ] }, { "cell_type": "code", "execution_count": null, "id": "3f80f939", "metadata": {}, "outputs": [], "source": [ "# set the random seed\n", "states, streams = sim_phtraj_from_times(opt_model,time_sim,seed=487)\n", "# generate a second path\n", "states, streams = sim_phtraj_from_times(opt_model,time_sim)\n", "states, streams" ] }, { "cell_type": "markdown", "id": "48d1414d", "metadata": {}, "source": [ "\n", "## 4.3 Recreating Data\n", "\n", "A final use-case is to use the simulation functions and the **times** of real data- this will let the user see if a given fit **model** properly describes the data (based on whatever statistics are relevant, eg. E and S).\n", "\n", "#### Simulating fit data based on burst times" ] }, { "cell_type": "code", "execution_count": null, "id": "372d98ee", "metadata": {}, "outputs": [], "source": [ "paths = []\n", "streams = []\n", "for tm in times:\n", " path_temp, stream_temp = sim_phtraj_from_times(opt_model,tm)\n", " paths.append(path_temp)\n", " streams.append(stream_temp)\n", "paths[0],streams[0]" ] }, { "cell_type": "markdown", "id": "3c408a44", "metadata": {}, "source": [ "\n", "# 5. Multiparameter H2MM\n", "\n", "\n", "## 5.1 Working with more than 2 streams\n", "\n", "The previous data used only 2 **streams** of data, so all the values in the `color` arrays were either 0 or 1. However, in many experiments, there will be more **streams**-- photons can be distinguished in a number of ways. *(think donor vs acceptor excitation in addition to emission, or 3 color FRET, or even polarization)*\n", "\n", "So how do we handle this? Well, pretty simple, each new stream gets a new number (an index to be precise) in the `color` array. So if there are $n$ streams, they will be indexed successively from $i=[0, 1, ..., n-1]$.\n", "\n", "So let's load a data set that has 3 **streams** (this is the same code as before, just a different file, and stored in the variables `color3` and `times3`)\n", "\n", "**For your data, it's your job to export it, if you're using *FRETBursts*, you should be able to copy-paste the `data_sort` function from the jupyter notebooks that go along with the mpH2MM Nature Communications paper.**\n", "\n", "#### Load 3 stream data" ] }, { "cell_type": "code", "execution_count": null, "id": "e3018c00", "metadata": {}, "outputs": [], "source": [ "times3 = []\n", "color3 = []\n", "i = 0\n", "# open the file\n", "with open('sample_data_3det.txt') as f:\n", " for line in f: # this reads each line one at a time\n", " if i % 2 == 0: # the datafiel is formate\n", " times3.append(np.array([int(x) for x in line.split()],dtype='Q'))\n", " else:\n", " color3.append(np.array([int(x) for x in line.split()],dtype='L'))\n", " i += 1" ] }, { "cell_type": "markdown", "id": "b00a0391", "metadata": {}, "source": [ "#### Let's inspect a color array to show that it has 3 **streams**" ] }, { "cell_type": "code", "execution_count": null, "id": "a03adabe", "metadata": {}, "outputs": [], "source": [ "color3[0]" ] }, { "cell_type": "markdown", "id": "47f4f273", "metadata": {}, "source": [ "See how instead of being a \"random\" set of 0's and 1's, it's now a \"random\" set of 0's, 1's and 2's? That means it's a 3 **stream** model.\n", "\n", "\n", "## 5.2 Setting up the **emission probability matrix**\n", "\n", "With 3 **streams** the **emission probability** matrix has 3 columns, instead of the 2 that it had with the previous 2 **stream** data. The **streams** and columns of a model must match, 3 **streams** 3 columns, 4 **streams** 4 columns. The **initial probability** and **transition probability** matrices are *unchanged*.\n", "\n", "To demonstrate that only **emission probability** matrix changes, we'll use the **initial probability** and **transition probability** matrices we defined all the way back in [section 1.4](#14) (so make sure you've run those cells, otherwise this next cell won't work)\n", "\n", "#### Define some new `.obs` arrays, and make 3 stream models" ] }, { "cell_type": "code", "execution_count": null, "id": "8021a589", "metadata": {}, "outputs": [], "source": [ "# define the obs arrays\n", "# note that these arrays all have 3 columns, meaning there are 3 streams in the data\n", "obs_1state3 = np.array([[0.3,0.3,0.4]]) \n", "obs_2state3 = np.array([[0.5,0.1,0.4],[0.1,0.5,0.4]])\n", "obs_3state3 = np.array([[0.4,0.4,0.2],[0.3,0.1,0.6],[0.2,0.4,0.4]])\n", "obs_4state3 = np.array([[0.4,0.4,0.2],[0.3,0.1,0.6],[0.2,0.4,0.4],[0.1,0.1,0.8]])\n", "\n", "# make the models\n", "init_model_1state3 = h2mm_model(prior_1state,trans_1state,obs_1state3)\n", "init_model_2state3 = h2mm_model(prior_2state,trans_2state,obs_2state3)\n", "init_model_3state3 = h2mm_model(prior_3state,trans_3state,obs_3state3)\n", "init_model_4state3 = h2mm_model(prior_4state,trans_4state,obs_4state3)" ] }, { "cell_type": "markdown", "id": "1ab2152b", "metadata": {}, "source": [ "#### Run the optimization\n", "\n", "Note that other than making sure the number of **streams** in the data matches the number of columns in the **emission probability** matrix, optimization and *Viterbi* work exactly the same." ] }, { "cell_type": "code", "execution_count": null, "id": "269b2dce", "metadata": {}, "outputs": [], "source": [ "opt_model_1state3 = EM_H2MM_C(init_model_1state3,color3,times3)\n", "opt_model_2state3 = EM_H2MM_C(init_model_2state3,color3,times3)\n", "opt_model_3state3 = EM_H2MM_C(init_model_3state3,color3,times3)\n", "opt_model_4state3 = EM_H2MM_C(init_model_4state3,color3,times3)" ] }, { "cell_type": "markdown", "id": "a06e0205", "metadata": {}, "source": [ "#### Run *Viterbi* to get the $ICL$" ] }, { "cell_type": "code", "execution_count": null, "id": "03ccb164", "metadata": {}, "outputs": [], "source": [ "path_1state3, scale_1state3, loglik_1state3, icl_1state3 = viterbi_path(opt_model_1state3,color3,times3)\n", "path_2state3, scale_2state3, loglik_2state3, icl_2state3 = viterbi_path(opt_model_2state3,color3,times3)\n", "path_3state3, scale_3state3, loglik_3state3, icl_3state3 = viterbi_path(opt_model_3state3,color3,times3)\n", "path_4state3, scale_4state3, loglik_4state3, icl_4state3 = viterbi_path(opt_model_1state3,color3,times3)" ] }, { "cell_type": "markdown", "id": "5e2bc9ac", "metadata": {}, "source": [ "#### Plot $BIC$" ] }, { "cell_type": "code", "execution_count": null, "id": "a9c1e7db", "metadata": {}, "outputs": [], "source": [ "opt_model_list3 = [opt_model_1state3,opt_model_2state3,opt_model_3state3,opt_model_4state3]\n", "bic3 = [model.bic for model in opt_model_list3]\n", "plt.scatter([1,2,3,4],bic3)\n", "plt.title(\"BIC\")" ] }, { "cell_type": "markdown", "id": "8c27ebaf", "metadata": {}, "source": [ "#### Plot $ICL$" ] }, { "cell_type": "code", "execution_count": null, "id": "2eac78ab", "metadata": {}, "outputs": [], "source": [ "icl3 = [icl_1state3,icl_2state3,icl_3state3,icl_4state3]\n", "plt.scatter([1,2,3,4],icl3)" ] }, { "cell_type": "markdown", "id": "92388604", "metadata": {}, "source": [ "So, it looks like the 3 state model is best for this set of data." ] }, { "cell_type": "markdown", "id": "fb246343", "metadata": {}, "source": [ "\n", "## 5.3 Checking a Model: `.ndet` and `.nstate` fields in `h2mm_model` objects\n", "\n", "Since `h2mm_model` objects have a number of states based on the size of all their matrices, as well as a number of **streams**, based on the number of columns in their **emission probability** matrix (`.obs`). \n", "\n", "We could check the number of states by checking the shape of `.prior` or `.trans`, even `.obs`, and we could check the number of detectors by checking the shape of `.obs`\n", "\n", "#### Check the states and streams of a model using the shapes of the `.trans` and `.obs` fields:" ] }, { "cell_type": "code", "execution_count": null, "id": "5ccb781d", "metadata": {}, "outputs": [], "source": [ "opt_model_2state3.trans.shape" ] }, { "cell_type": "markdown", "id": "4c0b9f8a", "metadata": {}, "source": [ "So this model has `.trans` of shape 2x2, so it is a 2 state model.\n", "\n", "We could find this with `opt_model_2state3.prior.shape` as well.\n", "\n", "If we look at the `.obs` field, we can get the number of **streams** as well:" ] }, { "cell_type": "code", "execution_count": null, "id": "2750c928", "metadata": {}, "outputs": [], "source": [ "opt_model_2state3.obs.shape" ] }, { "cell_type": "markdown", "id": "d8d7b86d", "metadata": {}, "source": [ "So the shape is `(2, 3)` which means that it is a 2 state model for 3 **streams**.\n", "\n", "#### There is a better way though:\n", "\n", "The `h2mm_model` object has two special fields: `.nstate` and `.ndet` for the number of **states** and **streams** respectively.\n", "#### Use `.nstate` and `.ndet` fields" ] }, { "cell_type": "code", "execution_count": null, "id": "e53e144d", "metadata": {}, "outputs": [], "source": [ "opt_model_2state3.nstate, opt_model_2state3.ndet" ] }, { "cell_type": "markdown", "id": "d1024380", "metadata": {}, "source": [ "These fields cannot be set, once a model is made, the number of **states** (`.nstate`) and **streams** is fixed (`.ndet`). These fields let us access these values in a streamlined way, which can be useful to make sure that a model has the appropriate number of streams for the data." ] }, { "cell_type": "markdown", "id": "0bef229b", "metadata": {}, "source": [ "\n", "# 6. Controlling `EM_H2MM_C`\n", "\n", "\n", "## 6.1 Convergence\n", "\n", "At it's core, **H2MM** is an optimization. Further, the mathematics and algorithm is formulated such that whenever the loglikelihood of a model is calculated, a new model is also generated, which is mathematically guaranteed to have a better *loglikelihood* than the previous model, unless the maximum has been reached. (The model is said then to have converged).\n", "\n", "Two (and sometimes related) issues arise:\n", "1. Floating point errors \n", " - The term \"maximum\" needs a little wiggle room, the models with differences in *loglikelihood* smaller than the floating point error should be considered as of equal quality.\n", "2. Long optimization\n", " - Coming to the \"maximum\" *loglikelihood* takes a very large number of iterations, and usually the differences both in values and *loglikelihoods* are minor.\n", " - This issue most often happens for large numbers of states, and overfit models.\n", "\n", "The solution comes by setting a maximum number of iterations the optimization will proceed, and setting a small, but non-zero threshold for the difference between *loglikelihoods* at which to consider sufficient to consider the optimization to have converged.\n", "\n", "By default, `EM_H2MM_C` will stop optimizing after 3600 iterations, and considers a difference in *loglikelihood* of $10^{-14}$ to be small enough to consider the optimization to have converged.\n", "\n", "These values can be overiden with the `max_iter` and `converged_min` keyword arguments (aka kwargs). \n", "\n", "*For Matlab users: keyword arguments are optional arguments with a default value. These are specified inside the parenthesis after a function, after all the normal arguments, and specified as `kwarg=value`*\n", "\n", "#### Changing `max_iter`" ] }, { "cell_type": "code", "execution_count": null, "id": "ab68193f", "metadata": {}, "outputs": [], "source": [ "small_opt_model_3state = EM_H2MM_C(init_model_3state, color, times, max_iter=500)" ] }, { "cell_type": "markdown", "id": "fce75c75", "metadata": {}, "source": [ "#### Changing `converged_min`" ] }, { "cell_type": "code", "execution_count": null, "id": "c7678eda", "metadata": { "scrolled": true }, "outputs": [], "source": [ "large_conv_model_3state = EM_H2MM_C(init_model_3state, color, times, converged_min=4e-6)" ] }, { "cell_type": "markdown", "id": "1a0b2517", "metadata": {}, "source": [ "There is a final keyword argument that can help control the length of optimization: `max_time`. It does almost as advertised, if the optimization goes over the time (which is in seconds), the optimization will stop at whatever model is currently being calculated.\n", "\n", "**HOWEVER** `max_time` *uses an **inaccurate clock**, which tends to run **fast** so it's use is not recommended*\n", "\n", "By default `max_time` is infinity, and so a model will never time out.\n", "\n", "\n", "### 6.1.1 Using `.optimize` instead\n", "\n", "The `.optimize` method takes the same keyword arguments as `EM_H2MM_C`, in the cell below, try using them to prove it to yourself." ] }, { "cell_type": "code", "execution_count": null, "id": "de7dbdfc", "metadata": {}, "outputs": [], "source": [ "# make a new model for yourself to use as a guinea pig\n", "play_model = h2mm_model(prior_3state,trans_3state,obs_3state3)\n", "# perform the optimization\n", "play_model.optimize(color3,times3,max_iter=800)" ] }, { "cell_type": "markdown", "id": "2b92ca80", "metadata": {}, "source": [ "\n", "## 6.2 Limits\n", "\n", "\n", "Sometimes we may want to restrain what values are possible in the H2MM model, for instance to keep transition rates within reasonable values, or because you know something about how the emmission probability matrix (`obs`) should behave.\n", "\n", "This is expected to happen most often (but not exclusively) when some experimental parameter is periodic, but not important to the data. For instance in $\\mu$sALEX experiments, the laser alternation period causes donor and acceptor excitation photons, which arrive in separate streams to alternate (and thus transition) perfectly periodically, yet that has no bearing on transitions between our system. Thus transition rates close to the rate of laser alternation are likely to be artifacts, and thus we want to exclude them. For $\\mu$sALEX experiments, we find this still is not enough, \n", "\n", "\n", "### 6.2.1 Our First Bounded Fit\n", "\n", "To define some bounds, we need to define the bounds, this is done using the `h2mm_limits` object, which we pass into the `EM_H2MM_C` function through the keyword argument `bounds`, and we also need to supply a method string to another keyword argument, `bounds_func`.\n", "\n", "Let's see an example, and then explain:\n", "#### Make an `h2mm_limits` object" ] }, { "cell_type": "code", "execution_count": null, "id": "57459f28", "metadata": {}, "outputs": [], "source": [ "alt_period = 4000 # a fake alternation period\n", "us_bounds = h2mm_limits(max_trans = 1/(alt_period))" ] }, { "cell_type": "markdown", "id": "df055d96", "metadata": {}, "source": [ "#### Run Optimization with bounds" ] }, { "cell_type": "code", "execution_count": null, "id": "ee6192c6", "metadata": {}, "outputs": [], "source": [ "# run optimization with bounds\n", "us_opt_model3 = EM_H2MM_C(init_model_3state3,color3,times3, bounds_func='revert',bounds=us_bounds)\n", "us_opt_model3" ] }, { "cell_type": "markdown", "id": "3f4f5441", "metadata": {}, "source": [ "So, what did we just do? The `h2mm_limits` object `us_bounds` prevents any value (off the diagonal) of the **transition probability** matrix (`.trans`) from ever being larger (i.e. faster transition rate) than `1/(4000)`. \n", "\n", "#### Bounds process\n", "When you use a bounds method, each iteration goes through the following steps:\n", "1. Calculate *loglikelihood* and new model\n", "2. Check if the **model** converged\n", "3. Analyze the **new model**, and correct if necessary\n", " 1. Check if any values are smaller or larger than a pre-set minimum or maximum\n", " 2. If values are out of bounds, apply correction, method defined by argument passed to `bounds_func`\n", "4. Repeat optimization (back to step 1)\n", "\n", "\n", "### 6.2.2 The `h2mm_limits` object\n", "\n", "`h2mm_limits` objects are designed specifically to be a flexible container for defining minimum and maximum values for elements of any and/or all of the `.prior`, `.trans`, and `.obs` arrays in an `h2mm_model` object during optimization.\n", "\n", "Notice that all arguments passed to `h2mm_limits()` are **keyword arguments**. Names are fairly obvious, with a `min_` or `max_` placed before the array name to define the minimum or maximum values respectively, so we have:\n", "- `min_prior` and `max_prior` to bound the `.prior` array\n", "- `min_trans` and `max_trans` to bound the `.trans` array\n", "- `min_obs` and `max_obs` to bound the `.obs` array\n", "\n", "If an argument is not specified, the default argument is either 0 or 1 for `min_` and `max_` which are the minimum and maximum possible in any **model**, so the bound imposes no practical limit to those values.\n", "\n", "These can be specified either as single floats (decimal number), or as arrays.\n", "\n", "Specifying as a float will mean that:\n", "1. All values will be bounded by the same minimum or maximum\n", " - For the `.trans` array, the values on the diagonal are ignored, as these will nearly always be close to 1\n", "2. The same `h2mm_limits`object can be used with `h2mm_model` objects with different **states** and **streams**\n", "\n", "Specifying as an array will mean that:\n", "1. Each element of the given array will be bounded separately, allowing greater precision\n", " - `.trans` diagonals are ignored\n", "2. The `h2mm_limits` object will only be able to be used with certain `h2mm_model` objects.\n", " - If any of the following `min_prior`, `max_prior`, `min_trans` or `max_trans` is specified as an array, then the `h2mm_limits` object must be used to optimize a `h2mm_model` object that has the same number of states\n", " - If `min_obs` or `max_obs` is specified, then the `h2mm_limits` object can only be used to optimize `h2mm_model` objects with the same number of **states** and **streams**\n", "\n", "In the previous example, we specified `max_trans` as a float. Therefore, it can be used to bound just about any model:\n", "\n", "#### Bounding several different models, with different numbers of streams and states:" ] }, { "cell_type": "code", "execution_count": null, "id": "5e59e641", "metadata": {}, "outputs": [], "source": [ "us_opt_model_1state = EM_H2MM_C(init_model_1state,color,times, bounds_func='minmax',bounds=us_bounds)\n", "us_opt_model_2state = EM_H2MM_C(init_model_2state,color,times, bounds_func='minmax',bounds=us_bounds)\n", "us_opt_model_3state = EM_H2MM_C(init_model_3state,color,times, bounds_func='minmax',bounds=us_bounds)\n", "us_opt_model_1state3 = EM_H2MM_C(init_model_1state3,color3,times3, bounds_func='minmax',bounds=us_bounds)\n", "us_opt_model_2state3 = EM_H2MM_C(init_model_2state3,color3,times3, bounds_func='minmax',bounds=us_bounds)\n", "us_opt_model_3state3 = EM_H2MM_C(init_model_3state3,color3,times3, bounds_func='minmax',bounds=us_bounds)" ] }, { "cell_type": "markdown", "id": "c7c06827", "metadata": {}, "source": [ "#### But, if we specify a different parameter with an array, then only matching models can be used:\n", "This time, let's work with the `min_obs` and `max_obs` arrays. This will restrict the individual states of the `.obs` matrix vary within certain ranges, each state to a different area. \n", "*This could be used for instance if we want to have only one state in a given area of the multi-dimensional space of the **emission probability** matrix, and all the other states in some other area of that range.*" ] }, { "cell_type": "code", "execution_count": null, "id": "c308f4df", "metadata": {}, "outputs": [], "source": [ "# define minimum and maximum arrays for obs, this will be a 4 state 3 detector model\n", "obs_min = np.array([[0.0,0.0,0.0],[0.1,0.0,0.3],[0.0,0.1,0.4],[0.0,0.0,0.7]])\n", "obs_max = np.array([[1.0,0.5,1.0],[0.6,0.3,0.7],[0.6,0.8,0.8],[1.0,1.0,1.0]])\n", "# make the limit object\n", "obs_limit = h2mm_limits(min_obs=obs_min, max_obs=obs_max)" ] }, { "cell_type": "markdown", "id": "c6edb5b5", "metadata": {}, "source": [ "#### Run the optimization" ] }, { "cell_type": "code", "execution_count": null, "id": "4e86b0cc", "metadata": {}, "outputs": [], "source": [ "obs_opt_model_4state3 = EM_H2MM_C(init_model_4state3,color3,times3,bounds=obs_limit,bounds_func='minmax')\n", "obs_opt_model_4state3" ] }, { "cell_type": "markdown", "id": "5f633f15", "metadata": {}, "source": [ "\n", "### 6.2.3 `bounds_func` Options\n", "\n", "There are several methods for \"correcting\" a new model that has \"out of bounds\" values. We can choose which one to use with the `bounds_func` keyword argument. These are all identified with strings.\n", "\n", "Our options are:\n", "- `'minmax'` when a value is out of bounds, it is set to the minimum or maximum depending on which was exceeded\n", "- `revert` resets the value to what it was in the current model (the one whose *loglikelihood* was just calculated\n", "- `revert_old` resets the value to the value in the model that preceded the current one in optimization\n", "\n", "How to choose between these three functions? It depends on what we want to limit/achieve. `'minmax'` will tend to just keep the values at their minimums/maximums, so it's great if we just don't want to optimize for really slow transition rates that are longer than our experiment, and thus never happen and therefore their values are more likely statistical noise. `'revert'` on the other hand will be more likely to be \"pushed off\" the path it was going and find a local maxima in the *loglikelihood*, away from artifacts produced by some experimental periodicity. But it also carries the greater risk of premature convergence. `'revert_old'` is basically just more extreme, even more likely to find a local maxima in *loglikelihood*, but also more likely to converge prematurely.\n", "\n", "#### Using an alternative `bounds_func` option\n", "*Try switching it between the different options*" ] }, { "cell_type": "code", "execution_count": null, "id": "dcc7a847", "metadata": {}, "outputs": [], "source": [ "usminmax_opt_model_2state3 = EM_H2MM_C(init_model_2state3,color3,times3, bounds_func='revert',bounds=us_bounds)\n", "usminmax_opt_model_2state3" ] }, { "cell_type": "markdown", "id": "75d0a784", "metadata": {}, "source": [ "Finally, be aware that `h2mm_limits` is very flexible, nearly any combination of min/max limits can be specified, and it can take mixing arrays and floats. The one limitation is that arrays must have matching **states** and **streams**. If conflicting arrays are given, an error will be raised, halting the Python script, and the conflict will need to be corrected. The `H2MM_C` API also checks for other potential conflicts in the value provided, i.e. a minimum cannot have a value greater than a maximum.\n", "\n", "#### A `h2mm_limits` object with multiple types of array specified mixing array and float definitions:" ] }, { "cell_type": "code", "execution_count": null, "id": "34044b34", "metadata": {}, "outputs": [], "source": [ "comb_bounds = h2mm_limits(max_trans=1/alt_period, min_obs=obs_min,max_obs=0.99999)\n", "comb_opt_model_2state3 = EM_H2MM_C(init_model_4state3,color3,times3, bounds_func='revert',bounds=comb_bounds)" ] }, { "cell_type": "markdown", "id": "1cb27f86", "metadata": {}, "source": [ "\n", "## 5.3 Auto-generate Initial Model with `factory_h2mm_model` function\n", "\n", "In [section 1.4](#wholeanalysis1), we had to define with all the arrays a new `h2mm_model` object for each number of **states** we wanted to optimize. This meant that if the system has 5 states, we'd have to define a 5th , and probably even a 6th model to be sure that 5 states is the best fit.\n", "\n", "It is better to define initial models based on prior knowledge of the system, but especially when we have a large number of states, or are trying to fit several different data sets, this can be rather laborious.\n", "\n", "The `H2MM_C` API contains a function to help us: `factory_h2mm_model`, which returns an initial `h2mm_model` object.\n", "\n", "The basic signature is: `factory_h2mm_model(nstate,ndet)` where `nstate` is the number of **states** in the output model, and `ndet` is the number of **streams** in the output model.\n", "\n", "\n", "### 6.3.1 Use `factory_h2mm_model` to make initial model\n", "*Try changing the arguments, and see how it changes the **streams** and **states** of the model.*" ] }, { "cell_type": "code", "execution_count": null, "id": "c237238e", "metadata": {}, "outputs": [], "source": [ "f_init_model_3state3 = factory_h2mm_model(3,3)\n", "f_init_model_3state3" ] }, { "cell_type": "markdown", "id": "8ebc2ad2", "metadata": {}, "source": [ "So let's re-write the code in [Section 1.4](#wholeanalysis1) using the `factory_h2mm_model`function:\n", "\n", "\n", "#### Basic Analysis with `factory_h2mm_model` to generate models" ] }, { "cell_type": "code", "execution_count": null, "id": "f494571c", "metadata": {}, "outputs": [], "source": [ "# load data (this is the same as before)\n", "times = []\n", "color = []\n", "i = 0\n", "# open the file\n", "with open('sample_data_2det.txt') as f:\n", " for line in f: # this reads each line one at a time\n", " if i % 2 == 0: # the datafiel is formate\n", " times.append(np.array([int(x) for x in line.split()],dtype='Q'))\n", " else:\n", " color.append(np.array([int(x) for x in line.split()],dtype='L'))\n", " i += 1\n", "\n", "\n", "# Define models using factory_h2mm_model\n", "finit_model_1state = factory_h2mm_model(1,2)\n", "finit_model_2state = factory_h2mm_model(2,2)\n", "finit_model_3state = factory_h2mm_model(3,2)\n", "finit_model_4state = factory_h2mm_model(4,2)\n", "\n", "\n", "# Optimize each model\n", "fopt_model_1state = EM_H2MM_C(finit_model_1state,color,times)\n", "fopt_model_2state = EM_H2MM_C(finit_model_2state,color,times)\n", "fopt_model_3state = EM_H2MM_C(finit_model_3state,color,times)\n", "fopt_model_4state = EM_H2MM_C(finit_model_4state,color,times)\n", "\n", "\n", "# calculate ICL\n", "path_1state, scale_1state, ll_1state, icl_1state = viterbi_path(fopt_model_1state,color,times)\n", "path_2state, scale_2state, ll_2state, icl_2state = viterbi_path(fopt_model_2state,color,times)\n", "path_3state, scale_3state, ll_3state, icl_3state = viterbi_path(fopt_model_3state,color,times)\n", "path_4state, scale_4state, ll_4state, icl_4state = viterbi_path(fopt_model_4state,color,times)" ] }, { "cell_type": "markdown", "id": "2a9c8bd7", "metadata": {}, "source": [ "\n", "### 6.3.2 Bounds and `factory_h2mm_model`\n", "\n", "The `factory_h2mm_model` function however can create problems when we are also bounding the values with `h2mm_limits` objects ([section 5.2](#52)). As it might generate a model with values that are out of the defined bounds. This will result in `EM_H2MM_C` potentially behaving oddly. This issue can be avoided by supplying the `h2mm_limits` object to `factory_h2mm_model` through the keyword argument `bounds`.\n", "\n", "#### Factory model with bounds" ] }, { "cell_type": "code", "execution_count": null, "id": "b0f409cb", "metadata": {}, "outputs": [], "source": [ "us_bounds = h2mm_limits(max_trans = 1/4000)\n", "binit_model_3state3 = factory_h2mm_model(3,3,bounds=us_bounds)\n", "binit_model_3state3" ] }, { "cell_type": "markdown", "id": "6d7427db", "metadata": {}, "source": [ "**The `h2mm_limits` object used for bounding should always be supplied to `factory_h2mm_model`**" ] }, { "cell_type": "markdown", "id": "5e87c349", "metadata": {}, "source": [ "\n", "### 6.3.3 Distribution kwargs in `factory_h2mm_model`\n", "\n", "There are a few additional keyword arguments to control how `factory_h2mm_model` distributes the values in the matrices of the output `h2mm_model` object.\n", "\n", "These are: \n", "- `prior_dist` valid values are `'even'` or `'random'`\n", "- `trans_dist` valid values are `'even'` or `'random'`\n", "- `obs_dist` valid values are `'even'` or `'random'`\n", "\n", "which set how values are distributed in the `.prior`, `.trans` and `.obs` arrays. There are only two possible inputs to these: `'even'` and `'random'` (both stings).\n", "\n", "When the value is set to `'even'` (the default), then values are evenly distributed inside the respective arrays.\n", "\n", "When `'random'` is chosen on the other hand, values will be assigned randomly in the respective arrays.\n", "\n", "`'even'` is usually a better choice, as the output is predictable, and an even distribution is generally a better starting point for optimization. `'random'` is useful for generating multiple initial **models** of the same **states** and **streams** to check if optimization finds the local or global maximum for the data. Further, occasionally when the `min_obs` or `max_obs` of the `bounds` are specified as arrays, setting `obs_dist` to `'random'` may be beneficial\n", "\n", "A final keyword argument, `trans_scale` is a float, and sets the off-diagonal values of the **transition probability** matrix (`.trans`). If `trans_dist = 'even'` and the default value does not conflict with the `bounds` input, then all off-diagonal values will be set to `trans_scale`, if `trans_dist = 'random'` then the values will be random, but centered around `trans_scale`.\n", "\n", "With all this information, it is now possible to re-write our analysis to limit the model transition rates for $\\mu$sALEX measurements, while also using the `factory_h2mm_model` function:\n", "\n", "\n", "### 6.3.4 Full analysis with bounds and `factory_h2mm_model`" ] }, { "cell_type": "code", "execution_count": null, "id": "d61dc500", "metadata": {}, "outputs": [], "source": [ "# load data (this is the same as before)\n", "times = []\n", "color = []\n", "i = 0\n", "# open the file\n", "with open('sample_data_3det.txt') as f:\n", " for line in f: # this reads each line one at a time\n", " if i % 2 == 0: # the datafiel is formate\n", " times3.append(np.array([int(x) for x in line.split()],dtype='Q'))\n", " else:\n", " color3.append(np.array([int(x) for x in line.split()],dtype='L'))\n", " i += 1\n", "\n", "\n", "# Define bounds object\n", "alt_period = 4000\n", "us_bounds = h2mm_limits(max_trans=1/alt_period)\n", "\n", "\n", "# Define models using factory_h2mm_model\n", "finit_model_1state3 = factory_h2mm_model(1,3,bounds=us_bounds)\n", "finit_model_2state3 = factory_h2mm_model(2,3,bounds=us_bounds)\n", "finit_model_3state3 = factory_h2mm_model(3,3,bounds=us_bounds)\n", "finit_model_4state3 = factory_h2mm_model(4,3,bounds=us_bounds)\n", "\n", "\n", "# Optimize each model\n", "fopt_model_1state3 = EM_H2MM_C(finit_model_1state3,color3,times3,bounds=us_bounds,bounds_func='revert')\n", "fopt_model_2state3 = EM_H2MM_C(finit_model_2state3,color3,times3,bounds=us_bounds,bounds_func='revert')\n", "fopt_model_3state3 = EM_H2MM_C(finit_model_3state3,color3,times3,bounds=us_bounds,bounds_func='revert')\n", "fopt_model_4state3 = EM_H2MM_C(finit_model_4state3,color3,times3,bounds=us_bounds,bounds_func='revert')\n", "\n", "\n", "# calculate ICL\n", "path_1state3, scale_1state3, ll_1state3, icl_1state3 = viterbi_path(fopt_model_1state3,color3,times3)\n", "path_2state3, scale_2state3, ll_2state3, icl_2state3 = viterbi_path(fopt_model_2state3,color3,times3)\n", "path_3state3, scale_3state3, ll_3state3, icl_3state3 = viterbi_path(fopt_model_3state3,color3,times3)\n", "path_4state3, scale_4state3, ll_4state3, icl_4state3 = viterbi_path(fopt_model_4state3,color3,times3)\n", "\n", "plt.scatter([1,2,3,4],[icl_1state3,icl_2state3,icl_3state3,icl_4state3])" ] }, { "cell_type": "markdown", "id": "899359d0", "metadata": {}, "source": [ "\n", "## 6.4 \"Continuous\" Analysis\n", "\n", "In the above example, all optimizations were run together, then all *Viterbi* analysis This is more clear and easier to read, however, it makes it harder to add more **states** to the analysis, so it is often more convenient to optimize, perform *Viterbi*, and then decide whether to optimize with one more **state**, or if we are confident that we have found the optimal number of **states** (usually by seeing that one or more models are over-fit).\n", "\n", "So bellow is an example of such code:\n", "\n", "\n", "### Optimize -> *Viterbi* workflow:" ] }, { "cell_type": "code", "execution_count": null, "id": "5c0e2db8", "metadata": {}, "outputs": [], "source": [ "# load data (this is the same as before)\n", "times = []\n", "color = []\n", "i = 0\n", "# open the file\n", "with open('sample_data_3det.txt') as f:\n", " for line in f: # this reads each line one at a time\n", " if i % 2 == 0: # the datafiel is formate\n", " times3.append(np.array([int(x) for x in line.split()],dtype='Q'))\n", " else:\n", " color3.append(np.array([int(x) for x in line.split()],dtype='L'))\n", " i += 1\n", "\n", "\n", "# Define bounds object\n", "alt_period = 4000\n", "us_bounds = h2mm_limits(max_trans=1/alt_period)\n", "\n", "\n", "# 1 state\n", "finit_model_1state3 = factory_h2mm_model(1,3,bounds=us_bounds)\n", "fopt_model_1state3 = EM_H2MM_C(finit_model_1state3,color3,times3,bounds=us_bounds,bounds_func='revert')\n", "path_1state3, scale_1state3, ll_1state3, icl_1state3 = viterbi_path(fopt_model_1state3,color3,times3)\n", "\n", "# 2 state\n", "finit_model_2state3 = factory_h2mm_model(2,3,bounds=us_bounds)\n", "fopt_model_2state3 = EM_H2MM_C(finit_model_2state3,color3,times3,bounds=us_bounds,bounds_func='revert')\n", "path_2state3, scale_2state3, ll_2state3, icl_2state3 = viterbi_path(fopt_model_2state3,color3,times3)\n", "\n", "# 3 state\n", "finit_model_3state3 = factory_h2mm_model(3,3,bounds=us_bounds)\n", "fopt_model_3state3 = EM_H2MM_C(finit_model_3state3,color3,times3,bounds=us_bounds,bounds_func='revert')\n", "path_3state3, scale_3state3, ll_3state3, icl_3state3 = viterbi_path(fopt_model_3state3,color3,times3)\n", "\n", "# 3 state\n", "finit_model_4state3 = factory_h2mm_model(4,3,bounds=us_bounds)\n", "fopt_model_4state3 = EM_H2MM_C(finit_model_4state3,color3,times3,bounds=us_bounds,bounds_func='revert')\n", "path_4state3, scale_4state3, ll_4state3, icl_4state3 = viterbi_path(fopt_model_4state3,color3,times3)\n", "\n", "plt.scatter([1,2,3,4],[icl_1state3,icl_2state3,icl_3state3,icl_4state3])" ] }, { "cell_type": "markdown", "id": "6e2e64af", "metadata": {}, "source": [ "\n", "## 6.5 While loop analysis\n", "\n", "The following code is a full analysis, using all the main techniques we have learned. It is \"smart\" in analyzing the fitting results so that it only calculates enough states to be *relatively* confident that we have found the ideal model.\n", "\n", "The basic algorithm is based on the following loop:\n", "1. Create model with $n$ states\n", "2. Optimize model\n", "3. Calcualte $ICL$ with *Viterbi*\n", "4. Compare all calculated models to test if the optimal model has been found\n", " 1. Should have at least 1 model with more states than optimal $ICL$\n", " 2. Should have minimum number of states calculated\n", " 3. Set a maximum number of states to calculate- as a way to prevent over-long analysis\n", "\n", "The following code also incorporates `h2mm_limits` and `factory_h2mm_model`.\n", "In practice, it is preferred to define 1 to 4 or 5 state models manually, and then switch to `factory_h2mm_model` to define larger models. Note that now results are stored in lists\n", "\n", "\n", "#### Full loop-based analysis" ] }, { "cell_type": "code", "execution_count": null, "id": "5d2c8b89", "metadata": {}, "outputs": [], "source": [ "# Load data file\n", "times3 = []\n", "color3 = []\n", "with open('sample_data_3det.txt') as f:\n", " for i, line in enumerate(f): # this reads each line one at a time\n", " if i % 2 == 0: # the datafiel is formate\n", " times3.append(np.array([int(x) for x in line.split()],dtype='Q'))\n", " else:\n", " color3.append(np.array([int(x) for x in line.split()],dtype='L'))\n", "\n", "# make bounds object\n", "alt_period = 4000\n", "us_bounds = h2mm_limits(max_trans=1/alt_period)\n", "\n", "# setup lists for storing results\n", "init_model = []\n", "opt_model = []\n", "path = []\n", "scale = []\n", "ll = []\n", "icl = []\n", "n = 1 # tracks number of states the model in the current iteration of the while loop should have\n", "min_states = 4 # minimum number of states thas must be calculated\n", "max_states = 8 # number of states beyond which to stop optimizing\n", "cont = True\n", "\n", "while cont:\n", " # make the n state model\n", " init_model_curr = factory_h2mm_model(n,3,bounds=us_bounds)\n", " init_model.append(init_model_curr) # append to list so we can inspect it later\n", " # optimize\n", " opt_model_curr = EM_H2MM_C(init_model_curr,color3,times3,bounds=us_bounds,bounds_func='revert')\n", " opt_model.append(opt_model_curr) # append to list so it can be examined at the end of the anlaysis\n", " # Viterbi\n", " path_curr, scale_curr, ll_curr, icl_curr = viterbi_path(opt_model_curr,color3,times3)\n", " # append results so they can be inspected later\n", " path.append(path_curr)\n", " scale.append(scale_curr)\n", " ll.append(ll_curr)\n", " icl.append(icl_curr)\n", " # Evaluate if we have sufficient knowledge to know if we have found the best fit model\n", " if n < min_states: # minimum number of models to optimize\n", " cont = True\n", " n += 1 # add 1 to n for next number of states to be calculated\n", " elif n > max_states: # don't optimize beyond some number of models\n", " cont = False\n", " elif np.min(icl) == icl_curr: # this says that the current icl is the minimum, therefore we do not have at least 1 overfit model\n", " cont = True\n", " n += 1 # add 1 to n for next number of states to be calculated\n", " else: # we have at least 1 overfit model, therefore we can stop optimizing\n", " cont = False\n", "\n", "# evaluating the results\n", "best_model_index = np.argmin(icl) # get index of lists of the \"best\" model, ie the one with the minimu ICL\n", "\n", "# plot the ICL graph\n", "plt.scatter(np.linspace(1,n,n),icl)\n", "opt_model[best_model_index] # final display the \"best model\"" ] }, { "cell_type": "markdown", "id": "b8d9d832", "metadata": {}, "source": [ "\n", "# 7. Customizing `EM_H2MM_C`\n", "\n", "The `EM_H2MM_C` function (and the `.optimize` method) has many options and ways to customize it's function. These are all specified through keyword arguments. So let's go through them.\n", "\n", "\n", "## 7.1 Changing parallelization: the `num_cores` kwarg\n", "\n", "The `EM_H2MM_C` algorithm is parallelized so it can use multiple cores on the computer. By default, the algorithm asks Python for the number of cores the machine has, and will limit itself to use only that many cores. But, for whatever reason, if we want to change this number, we have the keyword argument `num_cores`. \n", "\n", "Why would we want to change `num_cores`? Maybe we have a leave a core or two open for other applications because we have some other task running. Further, each system is a little different, and while the default should be close to optimum, we might be able to eek out a tiny bit better performance by changing `num_cores`\n", "\n", "Let's pretend we want to do a bunch of other things in other programs while `EM_H2MM_C` optimizes a model, so we'll use just one core, leaving the others free for our precious other programs. To do this, set `num_cores=1`.\n", "\n", "#### Optimization using only one core" ] }, { "cell_type": "code", "execution_count": null, "id": "5d989a88", "metadata": {}, "outputs": [], "source": [ "onecore_model_3state3 = EM_H2MM_C(init_model_3state3,color3,times3,num_cores=1)\n", "onecore_model_3state3" ] }, { "cell_type": "markdown", "id": "ec76c743", "metadata": {}, "source": [ "Try playing with the number of cores, and with different initial models, and see how it affects the time.\n", "\n", "*NOTE: changing the number of cores may cause the number of iterations before it converges to chnage, this is due to floating point errors. Check the actual results, and values of a 1 core vs n core optimization will differ by less than $10^{-10}$, the same can cause variation from optimization to optimization.*\n", "\n", "*In Jupyter notebooks, there are so-called \"magic\" commands, which start with % if applied to a single line, and %% to apply to an entire cell. The %%timeit magic can tell you how long it takes to execute a certain command, so it is included in the cell below to let you measure the time it takes for optimizations to run.* " ] }, { "cell_type": "code", "execution_count": null, "id": "0c1cab5f", "metadata": {}, "outputs": [], "source": [ "%%timeit\n", "onecore_model_3state3 = EM_H2MM_C(init_model_3state3,color3,times3,num_cores=1)" ] }, { "cell_type": "markdown", "id": "0523e719", "metadata": {}, "source": [ "`num_cores` is also a keyword argument in `viterbi_path`, `viterbi_sort`, and `H2MM_arr`. In all cases the underlying operation is the same.\n", "\n", "For the Python programmer: The default value for `num_cores = os.cpu_count()//2` The `//2` is because most systems are hyperthreaded, especially the higher powered systems researchers are likely to be working on. Therefore `os.cpu_count()` will return twice as many physical cores as there really are. If your system is not hyperthreaded, then you will want to manually specify `num_cores`. But I've found that depending on which OS/distro I'm using, how exactly it handles the parallelization differs. Therefore some experimenting with `num_cores` is likely to gain you a little speed up over the default.\n", "\n", "For the **C** programmer: All parallelization is handeled at the **C** level. For Linux, `H2MM_C` uses pthreads, while in Windows, the CreateThread function is used. The value passed to `num_cores` is the number of threads the algorithm will use. To Python programmers, understand that in **C**, threads behave much more like Python processes, in that they can run simultaneously on multiple cores." ] }, { "cell_type": "markdown", "id": "3a5781c1", "metadata": {}, "source": [ "\n", "## 7.2 Printing `EM_H2MM_C` iterations\n", "\n", "So, what if we want to track the progress of an optimization? If you've been paying close attention to the console/terminal window, you will have noticed that \n", "\n", "```\n", "Iteration 0 ...\n", "Iteration 1 ...\n", "Iteration 2 ...\n", "...\n", "```\n", "\n", "gets printed to this window whenever `EM_H2MM_C` is running. This helps keep track of the optimization, wthout cluttering up the Jupyter Notebook. \n", "\n", "This behavior can be modified however with the `print_func` keyword argument. With this we have a number of options, that let us print to the Jupyter Notebook if we want a record of the iterations, or even silence all printing of progress completely, bellow, see the list of acceptable inputs:\n", "\n", "\n", "### 7.2.1 Basic `print_func` options\n", "\n", "- `'console'` the default, prints to the **console** window, provides basic information on the progress of the optimization \n", "- `'all'` prints a representation of the whole model to be optimized in the next iteration (this will be very verbose This option is very verbose, use only if you want to see the exact steps of the optimization\n", "- `'diff'` Print the difference between the previous and current model logliklihoods and the current loglikelihood (this is the default if no `bounds` keyword argument is specified)\n", "- `'comp'` Print the old and current loglikelihoods\n", "- `'iter'` Prints only the iteration number- compact way to still track optimization progress in Jupyter Notebook\n", "- `None` Suppresses all printing of iteration information, to both **console** and **jupyter notebook**\n", "\n", "*If iterations show up in chunks when printing to a jupyter notebook, this is because the display to the Jupyter notebook is routed through a buffer, so what is printed to the notebook is often a few milliseconds behind, and multiple buffered messages get stored at once*\n", "\n", "### Below you have an example, try changing it to the different options and see how it changes the output" ] }, { "cell_type": "code", "execution_count": null, "id": "2ce1e8c6", "metadata": {}, "outputs": [], "source": [ "print_opt_model_2state3 = EM_H2MM_C(init_model_2state3,color3,times3,print_func='diff')" ] }, { "cell_type": "markdown", "id": "0afee4b2", "metadata": {}, "source": [ "\n", "### 7.2.2 Customized display using your own function\n", "\n", "There is one final type of argument that can be passed to `print_func`: a user-written function. When `EM_H2MM_C` detects that a function is passed to `print_func`, it will pass a specific set of arguments to that function:\n", "\n", "`print_func(niter, new_model, current_model, old_model, t_iter, t_total)`\n", "\n", "NOTE: It is not necessary, but recommended to keep these variable names in the function declaration.\n", "\n", "where:\n", "- `niter` is the number of iterations\n", "- `new_model` is a `h2mm_model` object that represents the next model to be optimized (**before** checking for out of bounds values) note that its `.loglik` will be irrelevant because it has not been calculated yet.\n", "- `current_model` is a `h2mm_model` object that represents the model whose `.loglik` was just calculated\n", "- `old_model` is a `h2mm_model` object that represents the model from the previous iteration.\n", "- `t_iter` is a float which is the time in seconds based on the **inaccurate C clock** that it took to calculate the latest iteration\n", "- `t_total` is a float which is the time in seconds based on the **inaccurate C clock** that the full optimization has taken\n", "\n", "We must write a function that accepts these arguments, and prints out whatever information we want to be displayed.\n", "\n", "NOTE: This function should be short, just whatever pre-processing is needed and then a `print` statement.\n", "\n", "#### Write a function to print the BIC of each iteration" ] }, { "cell_type": "code", "execution_count": null, "id": "c6a99962", "metadata": {}, "outputs": [], "source": [ "def bic_print(niter, new_model, current_model, old_model, t_iter, t_total):\n", " # the function signature needs to have the given order, but the names don't matter\n", " # however for readability, it's best just to keep the function signature exactly as is\n", " #print(f\"The current BIC is {current.bic}, and we've been going for a full {t_total} seconds\")\n", " print(f\"Iteration:{niter}, the BIC is {current_model.bic}, iteration time {t_iter} seconds\")\n", " # also, these functions will usually be really short, like this one, just a single print statement" ] }, { "cell_type": "markdown", "id": "e0cfce65", "metadata": {}, "source": [ "Now, to pass this function, we just need to pass it to the `print_func` keyword argument **without** parenthesis.\n", "\n", "*For Matlab users, when Python sees a function name without parentheses, it does not try to execute it, rather it tells Python to pass that function along, and it can be called later by the function it was passed to.*\n", "\n", "#### Pass the `bic_print` function to `print_func` keyword argument" ] }, { "cell_type": "code", "execution_count": null, "id": "3eb4d052", "metadata": { "scrolled": true }, "outputs": [], "source": [ "bic_print_opt_model_2state3 = EM_H2MM_C(init_model_2state3,color3,times3,print_func=bic_print)" ] }, { "cell_type": "markdown", "id": "4ddfacbf", "metadata": {}, "source": [ "\n", "## 7.3 User Defined Limits Function\n", "\n", "When I said that the `bounds` keyword argument had to be `None` or a `h2mm_limits` object, and `bounds_func` had to be one of `None`, `'minmax'`, `'revert'`, or `'revert_old'`, well there is one exception: a function can be passed to `bounds_func`.\n", "\n", "So let's first review how the optimization works:\n", "When you use a bounds method, each iteration goes through the following steps:\n", "1. Calculate *loglikelihood* and new model\n", "2. Check if model converged\n", "3. Analyze the new model, and correct if necessary\n", " 1. Check if any values are smaller or larger than a pre-set minimum or maximum\n", " 2. If values are out of bounds, apply correction, method defined by argument passed to `bounds_func`\n", "\n", "When `EM_H2MM_C` detects that a function has been passed to `bounds_func`, then it hands over step 3 to the function passed to `bounds_func`.\n", "\n", "This function must have the following signature:\n", "\n", "`bounds_func(new_model,current_model,old_model,bound)`\n", "\n", "where\n", "- `new_model` is a `h2mm_model` object representing the unoptimized model to be calculated in the next iteration, note that it' loglik has not been calculated, it will behave closer to an initial model.\n", "- `current_model` is a `h2mm_model` object representing the model of the current iteration, whose `.loglik` was just calculated\n", "- `old_model` is a `h2mm_model` object representing the model from the previous iteration\n", "- `bound` is the argument passed to the `bound` keyword argument\n", "\n", "and it must return a `h2mm_model` object.\n", "\n", "*As with passing a function to `print_func` it is not required, but recommended to keep the variable names*\n", "\n", "**WARNING:** be very careful in writing your own limits function- since just about any model can be passed back, if the adjusted model returned is not a logical adjustment, optimization with poorly written `bounds_func` functions will return meaningless results. \n", " \n", "You know the overworn adage- \"With great power comes great responsibility\" -- Uncle Ben\n", "\n", "#### Define a model value limiting function:" ] }, { "cell_type": "code", "execution_count": null, "id": "84fe3045", "metadata": {}, "outputs": [], "source": [ "def sample_bounds(new_model,current_model,old_model,bound):\n", " # like the print_func function, it's usually best to just keep the function signature the same\n", " # grab the obs matrix\n", " obs = new_model.obs\n", " # set first row of obs matrix to bound\n", " obs[0,:] = bound\n", " # change the obs matrix of the new model\n", " new_model.obs = obs\n", " # return the adjusted model\n", " return new_model" ] }, { "cell_type": "markdown", "id": "7096eb5b", "metadata": {}, "source": [ "#### Now run `EM_H2MM_C` with the `sample_bouds` function passed to `bounds_func`:" ] }, { "cell_type": "code", "execution_count": null, "id": "8eda7ff2", "metadata": {}, "outputs": [], "source": [ "bnd = np.array([0.09,0.01,0.9])\n", "smp_bnd_opt_model_3state3 = EM_H2MM_C(init_model_3state3,color3,times3,bounds_func=sample_bounds,bounds=bnd)" ] }, { "cell_type": "markdown", "id": "86c5eed3", "metadata": {}, "source": [ "\n", "# 8. Appendices\n", "\n", "\n", "\n", "## 8.1 `h2mm_model` fields\n", "\n", "Here is a list, and description of all the properties you can access.\n", "\n", "So you can look at the model through its different attributes (see [section 2.1](#21)):\n", "_(the following will produce `numpy` arrays as outputs)_\n", "- `.prior` the optimized initial probability\n", "- `.trans` the optimized transition probability matrix\n", "- `.obs` the optimized emission probability matrix\n", "\n", "You can see how many **states**, **streams** and *free parameters* are in the model with the following (see [section 5.3](#53)):\n", "- `.nstate` the number of states in the model, this is the 0 dimension of `.prior`, `.trans` and `.obs` arrays\n", "- `.ndet` the number of photon streams in the model, the 1 dimension of the `.obs` array\n", "- `.k` the number of free parameters- this is $q^2+q*(r-1)-1$ where $q$ equals `.nstate` and $r$ equals `.ndet`.\n", "\n", "There are a number of other properties accessible, only after a model has been optimized (see [section 1.3](#13)):\n", "- `.loglik` the calculated loglikelihood of the model, a `float`\n", "- `.bic` the Bayesian Information Criterion ($BIC$), a `float`\n", "- `.nphot` the number of photons in the data set on which the model was optimized, an `int`\n", "- `.niter` the number of iterations used to optimize the model, ie either the maximum number of iterations or the iteration at which the loglik improvement was smaller than the minimum difference specified in `converged_min`, and `int`.\n", "\n", "A few other properties also exist, and tell you about how the model was optimized or *loglikelihood* calculated:\n", "\n", "- The following return a boolean `True` or `False` \n", " - `.is_conv` returns True if the optimization reached maximum loglikelihood, and False otherwise\n", " - `.is_opt` returns True or false depending on whether the model has been optimized- False if the model is an \"initial\" model, or the *loglikelihood* calculation was done without optimization\n", " - `.is_calc`returns True if the loglik has been calculated against a set of data, and False if the model is an \"initial model\"\n", "- The remaining return an int or a string respectively \n", " - `.conv_code` returns an `int`, representing the status of the model, i.e. has it been optimized with **data** did it converge etc.\n", " - `.conv_str` returns a string the informs the user of if/how the particular object of the `h2mm_model` class finished it's optimization, this is basically a human readable interpretation of `.conv_code`\n", "\n", "The fields higlighted in green above were not introduced before in this tutorial. This is because they are either self explanatory, or because a use for them was not immediately obvious. The `H2MM_C` API was designed to be very full featured, so potentially useful fields were included in addition to fields \n", "\n", "### Rundown of fields we have not seen yet\n", "\n", "\n", "### 8.1.1 Statistical Fields\n", "\n", "Recall that the $BIC$ is defined as: $BIC = -2 LL + K \\ln(n)$. The $K$ and $n$ values are accessible directly via there one fields:\n", "- `.k` = $K$\n", "- `.nphot` = $n$\n", "\n", "At a core level, the `.bic` field is generated by calculating it from the `.loglik`, `.k` and `.nphot` fields. `.nphot` could also be useful as a way to test if two models were optimized against the same set of data. It is unlikely that two datasets will have exactly the same number of photons. **NOTE: you can be assured two models were optimized against different data if `.nphot` is different, however, the converse is not true. However, if `.nphot` is the same, it is highly likely (not guaranteed) that they were optimized against the same set of data.**\n", "\n", "\n", "#### Let's test this, if this is true, then we should be able to calculate the $BIC$ with `.loglik'`, `.k` and `.nphot` " ] }, { "cell_type": "code", "execution_count": null, "id": "456cb9dd", "metadata": {}, "outputs": [], "source": [ "-2* fopt_model_3state.loglik + fopt_model_3state.k * np.log(fopt_model_3state.nphot) == fopt_model_3state.bic" ] }, { "cell_type": "code", "execution_count": null, "id": "0b9e0cd3", "metadata": {}, "outputs": [], "source": [ "-2* fopt_model_3state.loglik + fopt_model_3state.k * np.log(fopt_model_3state.nphot)" ] }, { "cell_type": "markdown", "id": "0d11994b", "metadata": {}, "source": [ "\n", "### 8.1.2 `.niter`\n", "\n", "This field has no impact on any other value. It is just a record of how many iterations `EM_H2MM_C` took. Most optimizations converge to values that are within floating point error regardless of the initial model (as long as the initial models have the same number of **states**), but the number of iterations that they take will vary.\n", "\n", "#### Try different `h2mm_model` objects that have been created in this notebook, and check their `.niter` values" ] }, { "cell_type": "code", "execution_count": null, "id": "bf458823", "metadata": {}, "outputs": [], "source": [ "fopt_model_2state.niter, opt_model_4state.niter" ] }, { "cell_type": "markdown", "id": "4d780029", "metadata": {}, "source": [ "\n", "### 8.1.3 Convergence Info\n", "\n", "The remaining fields (replicated below) are all about identifying how a model was calculated, and if they have converged.\n", "We'll use these in [section 7.2.1](#721) to decide about optimizing more.\n", "- The following return a boolean `True` or `False` \n", " - `.is_conv` returns True or False depending on whether or not the model reached maximum *loglikelihood*, or if the maximum number iterations/time was reached, false for initial models, ones in the middle of optimization, or when their *loglikelihood* has been calculated without optimization\n", " - `.is_opt` returns True or false depending on whether the model has been optimized- False if the model is an \"initial\" model, or the *loglikelihood* calculation was done without optimization\n", " - `.is_calc`returns True if the loglik has been calculated against a set of data, and False if the model is an \"initial model\"\n", "- The remaining return an int or a string respectively \n", " - `.conv_code` returns an int, representing the status of the model, ie has it been optimized with **data** did it converge etc.\n", " - 0: initial models\n", " - 1: models in the middle of optimization whose *loglikelihood* has been calculated- this is generally only seen for models passed to the `print_func` function inside `EM_H2MM_C`\n", " - -1: models in the middle of optimization whose *loglikelihood* has NOT been calculated- this is generally only seen in the `new_model` argument passed to the `print_func` function inside `EM_H2MM_C`\n", " - 2: unoptimized models whose *loglikelihood* has been calculated, this occurs for models whose *loglikelihood* has been calculated using `H2MM_arr` or the `.evaluate` method.\n", " - 3: Models that converged based on `converged_min` bounds\n", " - 4: Models that were optimized the `max_iter` number of iterations\n", " - 5: models that were optimized for the `max_time` amount of time\n", " - 6: models where an error occurred in the subsequent iteration\n", " - `.conv_str` returns a string the informs the user of if/how the particular object of the `h2mm_model` class finished it's optimization, this is basically a human readable interpretation of `.conv_code`\n", "\n", "**NOTE**: The `.conv_code` of 6 is rare, and these models should be used carefully, as they are the last good model before an error occurred in the optimization, usually due to a floating point value being too close to zero, resulting in a divide by 0 error." ] }, { "cell_type": "markdown", "id": "c23b1c76", "metadata": {}, "source": [ "\n", "## 8.2 `EM_H2MM_C` and `.optimize` Arguments and Fields\n", "\n", "`EM_H2MM_C` has some optional keyword parameters, which fall into three categories:\n", "1. Processor settings (see [section 7.1](#71)): \n", " - `num_cores=os.cpu_count()//2` sets the number of processors to use in the calculation, changing this will not affect the result of the calculation, only the speed of computation. The defaul (`os/cpu_count()//2`) finds the number of cpu cores that python detects, and havles it, as python usually detects virtual cores, and thus it is assumed that the actual number of physical cores is half that. This however is only true for hyperthreaded systems, and thus if your system is otherwise, you may want to change this setting.\n", "2. Convergence criterions (see [section 6.1](#61)):\n", " - `max_iter=3600` the maximum number of iterations the optimization will continue unless another convergence criterion is met\n", " - `max_time=np.inf` the maximum amount of time the optimization will proceed, this is set to inf, so there will never be a timeout. It should be noted that this time is assesed by the C-level clock, which is often inaccurate, so while the time is supposedly in seconds, if this setting is set to a finite value, the optimization will likely terminate sooner than expected\n", " - `converged_min=1e-14` the minimum difference between models for an optimization to be considered converged, i.e. if the improvement in logliklihood is less than this, the optimization will terminate. This is the most common termination criteriod, as no futher optimization will proceed\n", " - `reset_niter=True` sets whether or not to set the number of iterations to 0 at the beginning of optimization. If set to `False`, then `max_iter` will trigger not after `max_iter` iterations have been perfomred in the current optimization, but rather it will start the counter at whatever number the input model has as `.niter`. Note also that if a data set with a different number of photons is passed, then regardless of whether `reset_niter` is set to `True` or `False`, the iteration counter will be reset.\n", "3. Printing options (see [section 7.2](#72)):\n", " - `print_func='console'` sets how the user is updated on the progress of each iteration.\n", "4. Bounding inputs: methods of limiting the values and adjusting the functionality of the H2MM optimization (see [section 6.2](#62) and [section 7.3](#73)):\n", " - `bounds_func=None` the function to run be run after each iteration of the optimzation, used to bound or otherwise change the next model in the iteration. Built-in functions identified with strings (disucssed later), if given a Python function, that function will be run (also discussed later)\n", " - `bounds=None` additional input to the bounds_func function (discussed later)\n", "\n", "\n", "### 8.2.1 The `reset_niter` argument\n", "\n", "It is possible to feed an already optimized model back into `EM_H2MM_C`. This might be done for instance if `EM_H2MM_C` reached the maximum number of iterations, but the model doesn't seem right, and so we want to let it optimize some more.\n", "\n", "By default, when `reset_niter` is `True`, `EM_H2MM_C` counts the number of iterations from 0, however, by setting `reset_niter` to `False` it will count iterations from the `.niter` value in the initial model passed to `EM_H2MM_C`. The same logic applies to the `.optimize` method of the `h2mm_model` object.\n", "\n", "`reset_niter=False` should be used when a model returned after reaching the maximum number of iterations is to be optimized further, and you want the output model to reflect the cumulative number of iterations. When this is done however, `max_iter` must be increased each time the model is resubmitted for further optimization, otherwise the optimization will stop after the first round.\n", "\n", "#### Keeping `.niter` after multiple rounds of optimization" ] }, { "cell_type": "code", "execution_count": null, "id": "6a486ce9", "metadata": {}, "outputs": [], "source": [ "opt_multi_model_3state3 = EM_H2MM_C(init_model_3state3,color3,times3,max_iter=100,reset_niter=False)\n", "opt_multi_model_3state3.niter" ] }, { "cell_type": "code", "execution_count": null, "id": "e4e2e5f6", "metadata": {}, "outputs": [], "source": [ "opt_multi_model_3state3 = EM_H2MM_C(opt_multi_model_3state3,color3,times3,max_iter=200,reset_niter=False)\n", "opt_multi_model_3state3.niter" ] }, { "cell_type": "markdown", "id": "a0a2c7d4", "metadata": {}, "source": [ "Note how the `.niter` field is smaller if we do the same, but with `reset_niter=True`" ] }, { "cell_type": "code", "execution_count": null, "id": "5d291eea", "metadata": {}, "outputs": [], "source": [ "opt_multi_model_3state3 = EM_H2MM_C(init_model_3state3,color3,times3,max_iter=110,reset_niter=True)\n", "opt_multi_model_3state3 = EM_H2MM_C(opt_multi_model_3state3,color3,times3,max_iter=100,reset_niter=True)\n", "opt_multi_model_3state3.niter" ] }, { "cell_type": "markdown", "id": "261e4dfa", "metadata": {}, "source": [ "So, you made it to the end, or did you skip everything?\n", "\n", "If you find any errors/odd behavior, please give your feedback.\n", "\n", "Thank you for using `H2MM_C`" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.4" } }, "nbformat": 4, "nbformat_minor": 5 }