Full reimplementation of H2MM, in C, with wrappers for Python and use from command line.

This code is substantially faster.

Author: harripd

.ipynb_checkpoints/Example-H2MM-Analysis-ImportH2MMcode-checkpoint.ipynb

@@ -0,0 +1,119 @@
+// File: C_H2MM.c
+// Author: Paul David Harris
+// Purpose: main wrapping functions to take burst data and submit to central H2MM algorithm
+// Date: 13 Feb 2021
+
+#ifdef linux
+#include <unistd.h>
+#elif _WIN32
+#include <windows.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "C_H2MM.h"
+
+#define TRUE 1
+#define FALSE 0
+
+// data pre-processor, takes the raw photon times and turns them into deltas
+// allocating new memory for the deltas, and attaching the delta to the appropriate phstream
+pwrs* get_deltas(unsigned long num_burst, unsigned long *burst_sizes, unsigned long long **burst_times, unsigned long **burst_det, phstream *b)
+{
+	size_t i, j;
+	size_t max_delta = 1;
+	size_t *deltas_temp;
+	size_t *inds_temp;
+	if ((burst_sizes == NULL) || (burst_times == NULL) || (burst_det == NULL) || (b == NULL))
+	{
+		printf("get_deltas(): One or more of the pointer arguments is NULL\n");
+		return NULL;
+	}
+	pwrs *powers = (pwrs*) calloc(1,sizeof(pwrs));
+	for ( i = 0; i < num_burst; i++)
+	{
+		deltas_temp = (size_t*) calloc(burst_sizes[i],sizeof(size_t));
+		inds_temp = (size_t*) calloc(burst_sizes[i],sizeof(size_t));
+		deltas_temp[0] = 0;
+		inds_temp[0] = burst_det[i][0];
+		for ( j = 1; j < burst_sizes[i]; j++)
+		{
+			if (burst_times[i][j] < burst_times[i][j-1])
+			{
+				printf("Photon %ld in burst %ld arrives out of order, please sanitize your data\n",j,i);
+				free(powers);
+				return NULL;
+			}
+			else if ( burst_times[i][j] == burst_times[i][j-1] ) deltas_temp[j] = 0;
+			else
+			{
+				deltas_temp[j] = (size_t) (burst_times[i][j] - burst_times[i][j-1]) - 1;
+				if ( deltas_temp[j] > max_delta) max_delta = deltas_temp[j];
+			}
+			inds_temp[j] = burst_det[i][j];
+		}
+		b[i].delta = deltas_temp;
+		b[i].det = inds_temp;
+		b[i].nphot = burst_sizes[i];
+		//~ printf("burst_sizes[%d]: %d\n",i,burst_sizes[i]);
+	}
+	powers->max_pow = max_delta + 1;
+	powers->pow_list = (size_t*) calloc(powers->max_pow,sizeof(size_t)); // allocate the memory for an unfilled powers list, index reduce will handle filling it in
+	return powers;
+}
+
+h2mm_mod* C_H2MM(unsigned long num_bursts, unsigned long *burst_sizes, unsigned long long **burst_times, unsigned long **burst_det, h2mm_mod *in_model, lm *limits)
+{
+	size_t i, j;
+	// alocate variables
+	phstream *b = (phstream*) calloc(num_bursts,sizeof(phstream)); // remember to free all b[n]->delta to prevent memory leak
+	// process burst arrays
+	size_t num_burst = (size_t) num_bursts;
+	//~ printf("Getting deltas\n");
+	pwrs *powers = get_deltas(num_burst,burst_sizes,burst_times,burst_det,b); // note: allocates the powers->pow_list array, remember to free powers->pow_list before free powers or b, also, the stride lengths and td/tv/tq are not assigned (should be 0 because of calloc)
+	//~ printf("Got powers\n");
+	in_model->nphot = 0;
+	for ( i = 0; i < num_bursts; i++) in_model->nphot += burst_sizes[i];
+	if (powers == NULL)
+	{
+		printf("You have an out of order photon\n");
+		return NULL;
+	}
+	for ( i = 0; i < num_bursts; i++)
+	{
+		for ( j = 0; j < b[i].nphot; j++)
+		{
+			if ( b[i].det[j] >= in_model->ndet)
+			{
+				printf("Your data has more photon streams than your h2mm model\n");
+				return in_model;
+			}
+		}
+	}
+	powers->sk = in_model->nstate;
+	powers->sj = powers->sk * in_model->nstate; // set strides, since these never change, we keep them the same 
+	powers->si = powers->sj * in_model->nstate;
+	powers->sT = powers->si * in_model->nstate;
+	powers->A = (double*) calloc(powers->sj * powers->max_pow,sizeof(double));
+	powers->Rho = (double*) calloc(powers->sT * powers->max_pow,sizeof(double));
+	// reduce indexes
+	//~ index_reduce_end(num_burst, b, powers);
+	// run main routine
+	//~ printf("Entering main routine\n");
+	h2mm_mod *out_model = compute_ess_dhmm(num_burst, b, powers, in_model, limits);
+	// copy out to out_model
+	//~ printf("\nPrior after copy is: \n");
+	//~ for ( i = 0; i < out->nstate; i++) printf("%f ",out_model->prior[i]);
+	//~ printf("\n");
+	// free memory
+	for ( i = 0; i < num_bursts; i++){
+		free(b[i].delta);
+		free(b[i].det);
+	}
+	free(b);
+	free(powers->pow_list);
+	free(powers->Rho);
+	free(powers->A);
+	free(powers);
+	return out_model;
+}

.ipynb_checkpoints/H2MM-DecayFitting-FRETBurstsModules-checkpoint.ipynb

@@ -0,0 +1,138 @@
+// File:  C_H2MM.h
+// 
+// fwd_back_photonbyphoton
+#ifdef _WIN32
+#include <windows.h>
+#endif
+typedef struct
+{
+	size_t nphot;
+	size_t *delta;
+	size_t *det;
+} phstream;
+
+struct temp
+{
+	long len_burst;
+	unsigned long long *times;
+	long *detectors;
+	struct temp *next;
+};
+typedef struct temp temps;
+
+typedef struct
+{
+	size_t nstate;
+	size_t ndet;
+	size_t nphot;
+	size_t niter;
+	size_t conv;
+	double *prior;
+	double *trans;
+	double *obs;
+	double loglik;
+} h2mm_mod;
+
+typedef struct
+{
+	phstream *phot;
+	size_t max_phot;
+	size_t *cur_burst;
+	size_t num_burst;
+	size_t sk; // number of states, indexing chosen to match
+	size_t sj; // square of the number of states
+	size_t si; // cube of the number of states
+	size_t sT; // fourth power of number of states
+	double *Rho;
+	double *A;
+	h2mm_mod *current;
+	h2mm_mod *new;
+#ifdef linux
+	pthread_mutex_t *h2mm_lock;
+#endif
+} fback_vals;
+
+// C_H2MM
+typedef struct
+{
+	size_t max_iter;
+	size_t num_cores;
+	double max_time;
+	double min_conv;
+} lm;
+
+
+// rho_calc
+typedef struct // pwrs is a structure that contains pointers to the A (transmat) and Rho arrays, plus information on dimensions and which values are to be calculated
+{
+	size_t max_pow;
+	size_t sT;
+	size_t si;
+	size_t sj; // note that this is also the stride for the delta t in the A array
+	size_t sk; // note that this is also the number of states
+	size_t tv;
+	size_t tq;
+	size_t td;
+	size_t *pow_list;
+	double *A;
+	double *Rho;
+} pwrs;
+
+// viterbi
+typedef struct
+{
+	size_t nphot;
+	size_t nstate;
+	double loglik;
+	size_t *path;
+	double *scale;
+} ph_path;
+
+
+typedef struct
+{
+	size_t si;
+	size_t sT;
+	double *A;
+	phstream *phot;
+	ph_path *path;
+	h2mm_mod *model;
+} vit_vals;
+
+
+// C_H2MM
+
+pwrs* get_deltas(unsigned long num_burst, unsigned long *burst_sizes, unsigned long long **burst_times, unsigned long **burst_det, phstream *b);
+
+h2mm_mod* C_H2MM(unsigned long num_bursts, unsigned long *burst_sizes, unsigned long long **burst_times, unsigned long **burst_det, h2mm_mod *in_model, lm *limits);
+
+// fwd_back_photonbyphoton
+#ifdef linux
+void* fwd_back_PhotonByPhoton(void* burst);
+#elif _WIN32
+DWORD WINAPI fwd_back_PhotonByPhoton(void* burst);
+#endif
+
+void h2mm_normalize(h2mm_mod *model_params);
+
+h2mm_mod* compute_ess_dhmm(size_t num_bursts, phstream *b, pwrs *powers, h2mm_mod *in, lm *limits);
+
+// rho_calc
+
+void* rhoulate(void *vals);
+
+void* rho_all(size_t nstate, double* transmat, pwrs *powers);
+
+// viterbi
+#ifdef linux
+void* viterbi_burst(void* in_vals);
+#elif _WIN32
+DWORD WINAPI viterbi_burst(void* in_vals);
+#endif
+int viterbi(unsigned long num_bursts, unsigned long *burst_sizes, unsigned long long **burst_times, unsigned long **burst_det, h2mm_mod *model, ph_path *path_array);
+
+// C_H2MM_interface
+
+temps* burst_read(char *fname, size_t *n);
+
+h2mm_mod* h2mm_read(char *fname);

.ipynb_checkpoints/H2MM-DecayFitting-MultiChannel-checkpoint.ipynb

@@ -0,0 +1,100 @@
+// File: C_H2MM_interface.c
+// Author: Paul David Harris
+// Purpose: Wrapper function for commandline interface with C_H2MM
+// Date: 03 Apr 2021
+
+#ifdef _WIN32
+#include <windows.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "C_H2MM.h"
+
+int main(int argc, char **argv)
+{
+	unsigned long long **times;
+	unsigned long **detectors;
+	long *len_bursts;
+	size_t i, j;
+	size_t num_burst = 0;
+	char *eptr;
+	temps *head;
+	//temps *next;                                                                                   // unreferenced??
+	temps *tmp;
+	phstream *b;
+	h2mm_mod *in_model;
+	lm *limits = (lm*) calloc(1,sizeof(lm));
+	limits->max_iter = 3600;
+	limits->max_time = INFINITY;
+	limits->min_conv = 1e-14;
+	limits->num_cores = 4;
+	int n = 0;
+	if (argc > 2)
+	{
+		//~ printf("Read bursts\n");	
+		head = burst_read(argv[1],&num_burst);
+		//~ printf("Read model\n");
+		in_model = h2mm_read(argv[2]);
+		if (argc > 3) limits->num_cores = (size_t) strtol(argv[3],&eptr,10);
+		if (argc > 4) limits->max_iter = (size_t) strtol(argv[4],&eptr,10);
+		if (argc > 5) limits->max_time = strtod(argv[4],&eptr);
+		if (argc > 6) limits->min_conv = strtod(argv[5],&eptr);
+	}
+	else
+	{
+		fprintf(stderr,"Too few arguments");
+		exit(EXIT_FAILURE);
+	}
+	// recast linked list of burst_read into
+	// first step is to allocate the necessary arrays of pointers and lengths 
+	times = (unsigned long long**)calloc(num_burst, sizeof(unsigned long long*));
+	detectors = (unsigned long**) calloc(num_burst, sizeof(unsigned long long*));
+	len_bursts = (long*) calloc(num_burst, sizeof(long));
+	// now follow the linked list, and put the pointers and lengths in the arrays
+	for( i = 0; i < num_burst; i++)
+	{
+		len_bursts[i] = head->len_burst;
+		times[i] = head->times;
+		detectors[i] = head->detectors;
+		tmp = head;
+		head = head->next;
+		free(tmp); // don't keep the linked list structure around
+		//~ printf("Freed burst %d\n",i);
+	}
+	//~ printf("There are %d bursts\n",num_burst);
+	b = (phstream*) calloc(num_burst,sizeof(phstream));
+	printf("main(): num_burst: %d\n", (int)num_burst);
+	printf("main(): in_model->ndet: %d\n", (int)in_model->ndet);
+	printf("Entering main algorithm\n");
+	h2mm_mod *out_model = C_H2MM(num_burst,len_bursts,times,detectors,in_model,limits);
+	if (out_model == NULL) printf("You have an out of order photon\n");
+	else if (out_model == in_model) printf("You have too many detectors in your data than allowed in your model\n");
+	else
+	{
+		if (out_model->conv == 1)
+			printf("Model converged after %ld iterations\n",out_model->niter);
+		else if (out_model->conv ==2)
+			printf("Maxiumum iterations reached\n");
+		else if (out_model->conv == 3)
+			printf("Maximum time reached\n");
+		else
+			printf("NAN error, returning last successful model\n");
+		printf("The final model is:\n");
+		printf("Prior:\n");
+		for ( i = 0; i < out_model->nstate; i++) printf("%f ",out_model->prior[i]);
+		printf("\nTrans:\n");
+		for ( i = 0; i < out_model->nstate; i++)
+		{
+			for ( j = 0; j < out_model->nstate; j++) printf("%e ",out_model->trans[out_model->nstate * i + j]);
+			printf("\n");
+		}
+		printf("Obs\n");
+		for ( i = 0; i < out_model->ndet; i++)
+		{
+			for ( j = 0; j < out_model->nstate; j++) printf("%f ",out_model->obs[out_model->nstate * i + j]);
+			printf("\n");
+		}
+	}
+}