Skip to content

Commit 4c18c73

Browse files
sagipesagipe
committed
Initial commit of files from post-doc
1 parent da3d87c commit 4c18c73

31 files changed

+1773
-9
lines changed

‎Cfind_first_smaller_sample_in_array2.c

+162
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
/*
2+
* Finds the first smaller sample in the given array
3+
*
4+
* INPUTS:
5+
* array : [double/single column vector] a sorted array
6+
* value : [double/single]
7+
*
8+
* OUTPUT:
9+
* index : [double] the
10+
*
11+
* Sagi Perel, 10/2012
12+
*/
13+
14+
#include "mex.h"
15+
#include "matrix.h"
16+
17+
void
18+
mexFunction(
19+
int num_output_args, // Number of left hand side (output) arguments
20+
mxArray *output_arg[], // Array of left hand side arguments
21+
int num_input_args, // Number of right hand side (input) arguments
22+
const mxArray *input_arg[]) // Array of right hand side arguments
23+
{
24+
double *array, value, *ptr;
25+
float *array_f;
26+
int array_length, this_array_length;
27+
mxClassID array_type;
28+
int debug=0;
29+
int i;
30+
int start_idx, end_idx, middle_idx;
31+
double epsilon = 0.00001;
32+
33+
//sanity check for input arguments
34+
if( num_input_args != 2) mexErrMsgTxt( "Cfind_first_smaller_sample_in_array2.c: wrong syntax: Cfind_first_smaller_sample_in_array2(array, value)\n");
35+
if( mxGetN(input_arg[0]) != 1 && mxGetM(input_arg[0]) != 1 ) mexErrMsgTxt( "Cfind_first_smaller_sample_in_array2.c: the array must be a vector\n");
36+
array_type = mxGetClassID(input_arg[0]);
37+
if( array_type != mxDOUBLE_CLASS && array_type != mxSINGLE_CLASS ) mexErrMsgTxt( "Cfind_first_smaller_sample_in_array2.c: the array must be a column vector of type double or float\n");
38+
39+
if( mxGetN(input_arg[1]) != 1) mexErrMsgTxt( "Cfind_first_smaller_sample_in_array2.c: the value must be a double\n");
40+
if( mxGetM(input_arg[1]) != 1) mexErrMsgTxt( "Cfind_first_smaller_sample_in_array2.c: the value must be a double\n");
41+
42+
array_length = (mxGetM(input_arg[0]) > mxGetN(input_arg[0])) ? mxGetM(input_arg[0]) : mxGetN(input_arg[0]);
43+
44+
//read input arguments:
45+
ptr = (double*) mxGetData( input_arg[1]);
46+
value = *ptr;
47+
if(debug)
48+
printf("array_length=%d,value=%f\n",array_length,value);
49+
50+
51+
//-read the signal and find the first
52+
switch(array_type)
53+
{
54+
case mxDOUBLE_CLASS:
55+
array = (double *) mxGetData( input_arg[0]);
56+
57+
//verify the array is sorted
58+
if(array[0] < array[array_length-1])
59+
mexErrMsgTxt( "Cfind_first_smaller_sample_in_array2.c: array should be sorted in decreasing order\n");
60+
61+
//check for edge condition- if the first element is smaller than the value
62+
if((array[0] - value)<0)
63+
{
64+
output_arg[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
65+
*mxGetPr(output_arg[0]) = 1;
66+
return;
67+
}
68+
69+
//shrink the array size by half every time
70+
start_idx = 0;
71+
end_idx = array_length-1;
72+
while(1)
73+
{
74+
75+
this_array_length = end_idx - start_idx + 1;
76+
if(debug)
77+
printf("start_idx=%d end_idx=%d (%d)",start_idx, end_idx, this_array_length);
78+
79+
if(this_array_length <= 20)
80+
{
81+
middle_idx = search_array_serial(array, value, start_idx, this_array_length);
82+
if(debug)
83+
printf("\nthis_array_length<=20: middle_idx=%d\n",middle_idx);
84+
if(middle_idx < 0)
85+
{
86+
output_arg[0] = mxCreateDoubleMatrix(0, 0, mxREAL);
87+
return;
88+
}else{
89+
90+
output_arg[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
91+
*mxGetPr(output_arg[0]) = middle_idx;
92+
return;
93+
}
94+
}
95+
96+
middle_idx = (this_array_length % 2 == 0 )? (start_idx+(this_array_length/2)) : (start_idx+((this_array_length-1)/2));
97+
if(debug)
98+
printf(" middle_idx=%d,val=%f\n",middle_idx,array[middle_idx]);
99+
//shrink the array size by half
100+
if(mxIsNaN(array[middle_idx]))
101+
{
102+
if(middle_idx == start_idx)
103+
{
104+
output_arg[0] = mxCreateDoubleMatrix(0, 0, mxREAL);
105+
return;
106+
}else{
107+
middle_idx = middle_idx -1;
108+
}
109+
}else{
110+
/*
111+
//if the middle_idx is exactly the value we are looking for- numerical issues might cause to miss it
112+
if(array[middle_idx]-epsilon <= value && array[middle_idx]+epsilon >= value)
113+
{
114+
output_arg[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
115+
*mxGetPr(output_arg[0]) = middle_idx+1;
116+
return;
117+
}
118+
*/
119+
120+
//else continue as usual
121+
if((array[middle_idx] - value) < 0)
122+
{
123+
end_idx = middle_idx;
124+
}else{
125+
start_idx = middle_idx;
126+
}
127+
}
128+
}
129+
130+
break;
131+
case mxSINGLE_CLASS:
132+
array_f = (float *) mxGetData( input_arg[0]);
133+
for(i=0; i<array_length; i++)
134+
{
135+
if((array_f[i] - value) < 0)
136+
{
137+
output_arg[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
138+
*mxGetPr(output_arg[0]) = (i+1);
139+
return;
140+
}
141+
}
142+
break;
143+
}
144+
145+
//create an empty output, to be returned if no such item exists
146+
output_arg[0] = mxCreateDoubleMatrix(0, 0, mxREAL);
147+
}
148+
149+
150+
int search_array_serial(double* array, double value, int start_idx, int array_length)
151+
{
152+
int i;
153+
for(i=start_idx; i<start_idx+array_length; i++)
154+
{
155+
if((array[i] - value) <= 0)
156+
{
157+
return (i+1);
158+
}
159+
}
160+
return -1;
161+
}
162+

‎Cfind_first_smaller_sample_in_array2.mexa64

12.2 KB
Binary file not shown.

‎Cfind_first_smaller_sample_in_array2.mexw64

8.5 KB
Binary file not shown.

‎compute_average_data_within_overlapping_bins.m

+59
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
function normalized_bin_averages = compute_average_data_within_overlapping_bins(time, data, bin_size, bin_start_times, normalizing_mean, normalizing_std)
2+
% function normalized_bin_averages = compute_average_data_within_overlapping_bins(time, data, bin_size, bin_start_times, normalizing_mean, normalizing_std)
3+
%
4+
% Computes the average value of the data matrix in bins - specified by bin_start_times
5+
% The value per bin will be: (bin_average - normalizing_mean) / normalizing_std
6+
% This function assumes that the bin sizes are all the same and specified by bin_size_sec
7+
%
8+
% INPUTS:
9+
% time: [array] with length N
10+
% data: [matrix][N x P] N data samples for P variables
11+
% bin_size : [double] in the same time units as time, specifying the bin width
12+
% bin_start_times : [array] with length L (must be L>2 bins), in the same time units as time, specifying the beginning time of every bin
13+
% [normalizing_mean]: [double] Default: 0. Mean to normalize the data by
14+
% [normalizing_std] : [double] Default: 1. STD to normalize the data by
15+
%
16+
% OUTPUTS:
17+
% normalized_bin_averages: [array] [L x P] with the normalized average data value in every bin.
18+
% NaN values are filled for bins with no data in them
19+
%
20+
% Sagi Perel, 10/2012
21+
22+
if(nargin < 4 || nargin > 6)
23+
error('compute_average_data_within_overlapping_bins: wrong number of input arguments provided');
24+
end
25+
if(~isvector(time) || ~isvector(bin_start_times))
26+
error('compute_average_data_within_overlapping_bins: input arguments must be vectors');
27+
end
28+
[N, P] = size(data);
29+
if(length(time) ~= N)
30+
error('compute_average_data_within_overlapping_bins: time and data must have the number of samples');
31+
end
32+
if(isvector(data))
33+
data = make_column_vector(data);
34+
end
35+
if(~isscalar(bin_size))
36+
error('compute_average_data_within_overlapping_bins: bin_size must be a double');
37+
end
38+
39+
num_bins = length(bin_start_times);
40+
if(num_bins < 2)
41+
error('compute_average_data_within_bins: bin_start_times must have at least 2 bins');
42+
end
43+
if(~exist('normalizing_mean','var'))
44+
normalizing_mean = 0;
45+
end
46+
if(~exist('normalizing_std','var'))
47+
normalizing_std = 1;
48+
end
49+
50+
normalized_bin_averages = nan(num_bins,P);
51+
for i=1:num_bins
52+
data_to_average = data( time > bin_start_times(i) & time <= (bin_start_times(i)+bin_size), :); % [M x P]
53+
if(~isempty(data_to_average))
54+
bin_average = nanmean(data_to_average);% [1 x P]
55+
normalized_bin_averages(i,:) = (bin_average - normalizing_mean) ./ normalizing_std;
56+
end
57+
end
58+
59+

‎compute_mutual_information.m

+119
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
function [I pxy px py edges_x edges_y nxy nx ny] = compute_mutual_information(x, y, x_partition_type, Lx, y_partition_type, Ly)
2+
% function [I pxy px py edges_x edges_y nxy nx ny] = compute_mutual_information(x, y, x_partition_type, Lx, y_partition_type, Ly)
3+
%
4+
% Computes the mutual information between time series x and time series y
5+
% using a uniform partition of the x-y plane to Lx*Ly bins.
6+
% The bins are of unequal sizes, so that the number of elements in every bin is roughly
7+
% equal (equal entropy).
8+
%
9+
% INPUTS:
10+
% x : [array] of length N
11+
% y : [array] of length N
12+
% x_partition_type: [string] 'uniform' : uniform marginal distibutions by using varying bin widths
13+
% 'equal_bins' : divide x to Lx bins of equal widths
14+
% 'unique' : count unique values (doesn't use Lx)
15+
% Lx: [int] number of bins, Lx <= N
16+
% y_partition_type: [string] 'uniform' : uniform marginal distibutions by using varying bin widths
17+
% 'equal_bins' : divide y to Ly bins of equal widths
18+
% 'unique' : count unique values (doesn't use Ly)
19+
% Ly: [int] number of bins, Ly <= N
20+
%
21+
% OUTPUTS:
22+
% I : [double] mutual information: I(X,Y) = sum_over_x( sum_over_y ( p(x,y)* log2( p(x,y)/ (p(x)*p(y)) ) ) )
23+
% pxy:
24+
% px :
25+
% py :
26+
% edges_x:
27+
% edges_y:
28+
%
29+
% Sagi Perel, 09/2012
30+
31+
if(nargin < 6 || nargin > 6)
32+
error('compute_mutual_information: wrong number if input arguments provided');
33+
end
34+
if(~isvector(x))
35+
error('compute_mutual_information: x must be a vector');
36+
end
37+
N = length(x);
38+
if(~isvector(y))
39+
error('compute_mutual_information: y must be a vector');
40+
elseif(length(y)~=N)
41+
error('compute_mutual_information: length(x)~=length(y)');
42+
end
43+
if(~exist('x_partition_type','var') || isempty(x_partition_type))
44+
x_partition_type = 'uniform';
45+
elseif(~ischar(x_partition_type))
46+
error('compute_mutual_information: x_partition_type must be a string');
47+
elseif(~any(strcmp(x_partition_type,{'uniform','equal_bins','unique'})))
48+
error('compute_mutual_information: unknown value for x_partition_type');
49+
end
50+
if(~strcmp(x_partition_type,'unique')) % x_partition_type,'unique' does not use Lx
51+
if(~isscalar(Lx))
52+
error('compute_mutual_information: Lx must be a scalar');
53+
elseif(Lx>N)
54+
error('compute_mutual_information: Lx > N');
55+
end
56+
end
57+
if(~exist('y_partition_type','var') || isempty(y_partition_type))
58+
y_partition_type = 'uniform';
59+
elseif(~ischar(y_partition_type))
60+
error('compute_mutual_information: y_partition_type must be a string');
61+
elseif(~any(strcmp(y_partition_type,{'uniform','equal_bins','unique'})))
62+
error('compute_mutual_information: unknown value for y_partition_type');
63+
end
64+
if(~strcmp(y_partition_type,'unique')) % y_partition_type,'unique' does not use Ly
65+
if(~isscalar(Ly))
66+
error('compute_mutual_information: Ly must be a scalar');
67+
elseif(Ly>N)
68+
error('compute_mutual_information: Ly > N');
69+
end
70+
end
71+
72+
if(all(isnan(x)) || all(isnan(y)))
73+
I = NaN;
74+
pxy = [];
75+
px = [];
76+
py = [];
77+
edges_x = [];
78+
edges_y = [];
79+
return;
80+
end
81+
82+
% the marginal distributions p(x), p(y) are approximated by histograms
83+
% which divide [min(x) max(x)] and [min(y) max(y)] according to x/y_partition_type
84+
switch(x_partition_type)
85+
case 'uniform'
86+
[edges_x nx] = sp_uniform_hist(x, Lx);
87+
case 'equal_bins'
88+
[edges_x nx] = sp_equal_bins_hist(x, Lx);
89+
case 'unique'
90+
[edges_x nx] = sp_unique_hist(x);
91+
otherwise
92+
error('compute_mutual_information: unknown value for x_partition_type');
93+
end
94+
switch(y_partition_type)
95+
case 'uniform'
96+
[edges_y ny] = sp_uniform_hist(y, Ly);
97+
case 'equal_bins'
98+
[edges_y ny] = sp_equal_bins_hist(y, Ly);
99+
case 'unique'
100+
[edges_y ny] = sp_unique_hist(y);
101+
otherwise
102+
error('compute_mutual_information: unknown value for y_partition_type');
103+
end
104+
105+
px = nx/sum(nx);
106+
py = ny/sum(ny);
107+
108+
% estimate the joint probability p(x,y)
109+
[nxy edges_x edges_y] = sp_joint_hist(x, y, x_partition_type, Lx, y_partition_type, Ly);
110+
pxy = nxy./sum(nxy(:)); % [Lx x Ly]
111+
112+
% compute the mutual information:
113+
% I(X,Y) = sum_over_x( sum_over_y ( p(x,y)* log2( p(x,y)/ (p(x)*p(y)) ) ) )
114+
px_py = make_column_vector(px)*make_row_vector(py); % [Lx x Ly]
115+
% don't include zero values
116+
lidx = px_py(:)>0 & pxy(:)>0; % logical vector with length [Lx * Ly]
117+
% compute the sum over all x
118+
I = sum( pxy(lidx) .* ( log2(pxy(lidx)) - log2(px_py(lidx)) ) );
119+

‎sp_array_to_matrix.m

+51
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
function matrix = sp_array_to_matrix(array, indices_matrix)
2+
% function matrix = sp_array_to_matrix(array, indices_matrix)
3+
%
4+
% Maps values from array to a matrix, based on the provided indices
5+
%
6+
% INPUTS:
7+
% array : [array] with data values
8+
% indices_matrix: [matrix] with indices of data in array. NaN entries stay NaN in the output matrix.
9+
%
10+
% OUTPUTS:
11+
% matrix: [matrix] same size as indices_matrix, with data from array mapped to the right places
12+
%
13+
% EXAMPLES:
14+
% sp_array_to_matrix(1:9, [1 2 NaN 3 4; 6 7 8 9 5])
15+
% ans =
16+
% 1 2 NaN 3 4
17+
% 6 7 8 9 5
18+
%
19+
% sp_array_to_matrix(1:9, [1 1 NaN 1 1; 6 7 8 9 5])
20+
% ans =
21+
% 1 1 NaN 1 1
22+
% 6 7 8 9 5
23+
% Sagi Perel, 12/2012
24+
25+
if(nargin ~= 2)
26+
error('sp_array_to_matrix: wrong number of input arguments provided');
27+
end
28+
if(~isvector(array))
29+
error('sp_array_to_matrix: array should be a vector');
30+
else
31+
num_values = length(array);
32+
end
33+
if(~ismatrix(indices_matrix))
34+
error('sp_array_to_matrix: indices_matrix should be a matrix');
35+
else
36+
linearized_indices = indices_matrix(:);
37+
if(nanmin(linearized_indices)<1)
38+
error('sp_array_to_matrix: indices_matrix should contain positive integers');
39+
elseif(nanmax(linearized_indices)>num_values)
40+
error('sp_array_to_matrix: indices_matrix contains an index too large for array');
41+
elseif(sum(~isnan(linearized_indices)) ~= num_values)
42+
warning('sp_array_to_matrix: some values in array are not used');
43+
end
44+
end
45+
% linearized_indices has some size, with possible NaN values
46+
% so find indices of elements in linearized_indices which are not NaN
47+
not_NaN_lidxs = ~isnan(linearized_indices);
48+
array_mapped_to_linearized_indices = nan(size(linearized_indices));
49+
array_mapped_to_linearized_indices(not_NaN_lidxs) = array(linearized_indices(not_NaN_lidxs));
50+
matrix = reshape(array_mapped_to_linearized_indices, size(indices_matrix));
51+

0 commit comments

Comments
 (0)