text:
function partitions = cosmo_oddeven_partitioner(ds, type)
% generates an odd-even partition scheme
%
% partitions=cosmo_oddeven_partitioner(chunks,[type])
%
% Input
%    ds              dataset struct with field .ds.chunks.
%    type            One of:
%                    - 'full': two partitions are returned, training on odd
%                       unique chunks and testing on even unique chunks,
%                       and vice versa (default). A typical use case is
%                       classification analysis (using
%                       cosmo_crossvalidation_measure)
%                    - 'half': a single partition is returned, training on
%                       odd unique chunks and testing on even unique chunks
%                       only. A typical use case is split-half
%                       correlation analysis (using
%                       cosmo_correlation_measure) because correlations are
%                       symmetric (i.e. corr(a,b)==corr(b,a) for column
%                       vectors a and b)
%
% Output:
%    partitions      A struct with fields .train_indices and .test_indices.
%                    Each of these is an Nx1 cell (N=1 when type='half',
%                    N=2 when type='full'.).
%                    .train_indices{k} and .test_indices{k} contain the
%                    sample indices for the sets of unique chunks
%                    alternatingly
%
% Example:
%     ds=struct();
%     ds.samples=NaN(8,99); % will be ignored by this function
%     ds.sa.targets=[8 9 8 9 8 9 8 9]';
%     ds.sa.chunks=[1 1 2 2 6 7 7 6]';
%     p=cosmo_oddeven_partitioner(ds);
%     % note that chunks=6 ends up in the odd chunks and chunks=7 in the
%     % even chunks, as 6 [7] is the third [fourth] unique value of chunks.
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 1    [ 3
%     %||       2      4
%     %||       5      6
%     %||       8 ]    7 ] }
%     %|| .test_indices
%     %||   { [ 3    [ 1
%     %||       4      2
%     %||       6      5
%     %||       7 ]    8 ] }
%     %||
%     %
%     % only half-partition (for correlation-based analysis)
%     p=cosmo_oddeven_partitioner(ds,'half');
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 1
%     %||       2
%     %||       5
%     %||       8 ] }
%     %|| .test_indices
%     %||   { [ 3
%     %||       4
%     %||       6
%     %||       7 ] }
%
% Notes:
%   - typically, the 'half' option can be used with
%     cosmo_correlation_measure, because correlations are symmetric; the
%     'full' option can be used with cosmo_crossvalidation_measure.
%   - this function returns partitions based on the sorted unique values
%     in chunks, not the chunk values themselves. For example, if the
%     sorted unique values of .sa.chunks are [2,4,5,8], then the values
%     2 and 5 are at the 'odd' position (1 and 3), and 4 and 8 are at the
%     'even' position.
%
% See also: cosmo_nchoosek_partitioner
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    chunks=get_chunks(ds);

    if nargin<2 || isempty(type)
        type='full';
    end

    switch type
        case 'full'
            do_half_partition=false;
        case 'half'
            do_half_partition=true;
        otherwise
            error('illegal type: must be ''full'' or ''half''');
    end

    indices=cosmo_index_unique(chunks);
    nparts=numel(indices);
    if nparts<2
        error('Need >=2 chunks, found %d', numel(indices));
    end

    % there are two partitions, unless do_half_partition in which case
    % there is one

    if do_half_partition
        npartitions=1;
    else
        npartitions=2;
    end

    % allocate space for output
    train_indices=cell(1,npartitions);
    test_indices=cell(1,npartitions);

    % Make partitions using even and odd chunks
    % >@@>

    % find the indices of even and odd chunks
    odd_indices=cat(1,indices{1:2:end});
    even_indices=cat(1,indices{2:2:end});

    % set the train and test indices
    train_indices{1}=odd_indices;
    test_indices{1}=even_indices;

    if npartitions==2
        train_indices{2}=even_indices;
        test_indices{2}=odd_indices;
    end

    % <@@<

    partitions.train_indices=train_indices;
    partitions.test_indices=test_indices;

    % final check to make sure all is kosher (including balanced-ness)
    cosmo_check_partitions(partitions,ds);

function ds=get_chunks(ds)
    if isnumeric(ds) && isvector(ds)
        % direct numeric
        return
    elseif isstruct(ds)
        if cosmo_isfield(ds,'sa.chunks')
            ds=ds.sa.chunks;
            return
        end
    end

    error(['illegal input: expected dataset struct with field '...
                '.sa.chunks, or numeric vector']);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function output = base64decode(input)
%BASE64DECODE Decode Base64 string to a byte array.
%
%    output = base64decode(input)
%
% The function takes a Base64 string INPUT and returns a uint8 array
% OUTPUT. JAVA must be running to use this function. The result is always
% given as a 1-by-N array, and doesn't retrieve the original dimensions.
%
% See also base64encode

error(javachk('jvm'));
if ischar(input), input = uint8(input); end

output = typecast(org.apache.commons.codec.binary.Base64.decodeBase64(input), 'uint8')';

end

author: happynear

repo: FaceVerification


text:
% NANMAX provides a replacement for MATLAB's nanmax that is almost
% compatible.

function [varargout] = nanmax(varargin)

% the standard MATLAB max function ignores NaNs from at least version 7.0.4 onwards
[varargout{1:nargout}] = max(varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function a=IGHM(b)
%IGHM Single-level inverse discrete 2-D multiwavelet transform.
%   IGHM performs a single-level 2-D multiwavelet reconstruction
%   using GHM multiwavelet with four multi-filters
%
%   Y = IGHM(X) computes the original matrix from the approximation
%   coefficients matrix LL and details coefficients matrices
%   LH, HL, HH, obtained by a multiwavelet decomposition of the
%   original input matrix and puts the result in Y.
%
%   The size of Y is half that of X which should be
%   a square matrix of size NxN where N is power of 2 since
%   X is de-vectorized by a repeated row postprocessing.
%
%   X should be arranged as [LL,LH;HL,HH].
%
%   See also GHM, GHMAP, IGHMAP DD2, IDD2, DWTMODE, WAVEDEC2, WAVEINFO.

%   Auth: Dr. Bessam Z. Hassan
%   Last Revision: 27-Feb-2004.
%   Copyright 1995-2002 The MathWorks, Inc.
% $Revision: 1.0 $

% initialize the coefficients
H0=[3/(5*sqrt(2)),4/5;-1/20,-3/(10*sqrt(2))];
H1=[3/(5*sqrt(2)),0;9/20,1/sqrt(2)];
H2=[0,0;9/20,-3/(10*sqrt(2))];
H3=[0,0;-1/20,0];
G0=[-1/20,-3/(10*sqrt(2));1/(10*sqrt(2)),3/10];
G1=[9/20,-1/sqrt(2);-9/(10*sqrt(2)),0];
G2=[9/20,-3/(10*sqrt(2));9/(10*sqrt(2)),-3/10];
G3=[-1/20,0;-1/(10*sqrt(2)),0];
% check the input
if nargin == 0,
	error('Not enough input arguments.');
end
if isempty(b)
   a = [];
   return
end
[N,M]=size(b);N=ceil(N/2);M=ceil(M/2);
if N~=2^round(log(N)/log(2))
    error('size of the input should be power of 2');
end
if M~=N
    error('the input matrix must be square');
end
% construct the W matrix
w=[H0,H1,H2,H3;G0,G1,G2,G3];
for i=1:N/2-1
    W(4*(i-1)+1:4*i,4*i-3:4*i+4)=w;
end
W=[W;[[H2,H3;G2,G3],zeros(4,2*N-8),[H0,H1;G0,G1]]];
p=[];X=[];
% coefficient shuffling
b1=b(1:N,1:N);b2=b(1:N,N+1:2*N);b3=b(N+1:2*N,1:N);b4=b(N+1:2*N,N+1:2*N);
b1=b1';b2=b2';b3=b3';b4=b4';
T(1:2:N,:)=b1(1:N/2,:);T(2:2:N,:)=b1(N/2+1:N,:);T=T';b1(1:2:N,:)=T(1:N/2,:);b1(2:2:N,:)=T(N/2+1:N,:);
T(1:2:N,:)=b2(1:N/2,:);T(2:2:N,:)=b2(N/2+1:N,:);T=T';b2(1:2:N,:)=T(1:N/2,:);b2(2:2:N,:)=T(N/2+1:N,:);
T(1:2:N,:)=b3(1:N/2,:);T(2:2:N,:)=b3(N/2+1:N,:);T=T';b3(1:2:N,:)=T(1:N/2,:);b3(2:2:N,:)=T(N/2+1:N,:);
T(1:2:N,:)=b4(1:N/2,:);T(2:2:N,:)=b4(N/2+1:N,:);T=T';b4(1:2:N,:)=T(1:N/2,:);b4(2:2:N,:)=T(N/2+1:N,:);
b=[b1,b2;b3,b4];
p=b';
% column vector shuffling
ii=0:4:2*N-1;jj=sort([ii+1,ii+2]);kk=sort([ii+3,ii+4]);
z(jj,:)=p(1:N,:);z(kk,:)=p(N+1:2*N,:);
% column inverse transform
X=W'*z;
% downsampling columns (repeated row postprocessing)
aa=X(1:2:2*N,:);
p=aa';
% row vector shuffling
ii=0:4:2*N-1;jj=sort([ii+1,ii+2]);kk=sort([ii+3,ii+4]);z=[];
z(jj,:)=p(1:N,:);z(kk,:)=p(N+1:2*N,:);
% row inverse transform
X=W'*z;
%downsampling rows (repeated row postprocessing)
a=X(1:2:2*N,:);

author: Sable

repo: mcbench-benchmarks


text:
%% Machine Learning Online Class - Exercise 3 | Part 1: One-vs-all

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     lrCostFunction.m (logistic regression cost function)
%     oneVsAll.m
%     predictOneVsAll.m
%     predict.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Setup the parameters you will use for this part of the exercise
input_layer_size  = 400;  % 20x20 Input Images of Digits
num_labels = 10;          % 10 labels, from 1 to 10
                          % (note that we have mapped "0" to label 10)

%% =========== Part 1: Loading and Visualizing Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  You will be working with a dataset that contains handwritten digits.
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

load('ex3data1.mat'); % training data stored in arrays X, y
m = size(X, 1);

% Randomly select 100 data points to display
rand_indices = randperm(m);
sel = X(rand_indices(1:100), :);

% to display example digits, uncomment the below
%displayData(sel);
%fprintf('Program paused. Press enter to continue.\n');
%pause;

%% ============ Part 2: Vectorize Logistic Regression ============
%  In this part of the exercise, you will reuse your logistic regression
%  code from the last exercise. You task here is to make sure that your
%  regularized logistic regression implementation is vectorized. After
%  that, you will implement one-vs-all classification for the handwritten
%  digit dataset.
%

fprintf('\nTraining One-vs-All Logistic Regression...\n')

lambda = 0.1;
[all_theta] = oneVsAll(X, y, num_labels, lambda);

%fprintf('Program paused. Press enter to continue.\n');
%pause;

%% ================ Part 3: Predict for One-Vs-All ================
%  After ...
pred = predictOneVsAll(all_theta, X);

fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);

author: worldveil

repo: coursera-ml


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

clear; clc;
%% Parameters
S = 100;    % Spot prices
T = 5;      % maturities
r = 0;      % discount factors
d = 0;      % dividends

titS = 'MC NIG - Asset Paths';

legend_base = 'Base scenario';

alpha = 10; % parameter
beta = 0;   % parameter
delta = 1;  % parameter
mu = 0;     % parameter

%% Simulation parameters
NTime = 120; NSim = 1; NBatches = 1;
K = ones(NTime+1,1);
K = T*cumsum(K)/(NTime+1);
K1 = K(2:end);

a_par = T/NTime;                         % parameter of IG distribution
b_par = delta*sqrt(alpha^2-beta^2);   % parameter of IG distribution
theta = a_par/b_par;
chi = a_par^2;

rstream = RandStream('mt19937ar','Seed',12345);
rstreamstate = rstream.State;

PathS = MC_NIG(S,r,d,T,alpha,beta,delta,NTime,NSim,NBatches);
%% Changing lambda
alpha_low = 5;
alpha_high = 15;

    rstream.State = rstreamstate;
    PathS_low = MC_NIG(S,r,d,T,alpha_low,beta,delta,NTime,NSim,NBatches);
    rstream.State = rstreamstate;
    PathS_high = MC_NIG(S,r,d,T,alpha_high,beta,delta,NTime,NSim,NBatches);

legend_low = 'Changing \alpha low';
legend_high = 'Changing \alpha high';

createfigure_path(K,PathS_low, PathS, PathS_high, titS, legend_low, legend_base, legend_high);

% calculate the returns
ReturnsS1_low = calcreturns(PathS_low,0);
ReturnsS1 = calcreturns(PathS,0);
ReturnsS1_high = calcreturns(PathS_high,0);

createfigure_returns(K1,ReturnsS1_low, ReturnsS1, ReturnsS1_high, titS, legend_low, legend_base, legend_high);
%% Changing kappa
beta_low = -2;
beta_high = 2;

    rstream.State = rstreamstate;
    PathS_low = MC_NIG(S,r,d,T,alpha,beta_low,delta,NTime,NSim,NBatches);
    rstream.State = rstreamstate;
    PathS_high = MC_NIG(S,r,d,T,alpha,beta_high,delta,NTime,NSim,NBatches);

legend_low = 'Changing \beta low';
legend_high = 'Changing \beta high';

createfigure_path(K,PathS_low, PathS, PathS_high, titS, legend_low, legend_base, legend_high);

% calculate the returns
ReturnsS1_low = calcreturns(PathS_low,0);
ReturnsS1 = calcreturns(PathS,0);
ReturnsS1_high = calcreturns(PathS_high,0);

createfigure_returns(K1,ReturnsS1_low, ReturnsS1, ReturnsS1_high, titS, legend_low, legend_base, legend_high);
%% Changing eta
delta_low = 0.25;
delta_high = 1.75;

    rstream.State = rstreamstate;
    PathS_low = MC_NIG(S,r,d,T,alpha,beta,delta_low,NTime,NSim,NBatches);
    rstream.State = rstreamstate;
    PathS_high = MC_NIG(S,r,d,T,alpha,beta,delta_high,NTime,NSim,NBatches);

legend_low = 'Changing \delta low';
legend_high = 'Changing \delta high';

createfigure_path(K,PathS_low, PathS, PathS_high, titS, legend_low, legend_base, legend_high);

% calculate the returns
ReturnsS1_low = calcreturns(PathS_low,0);
ReturnsS1 = calcreturns(PathS,0);
ReturnsS1_high = calcreturns(PathS_high,0);

createfigure_returns(K1,ReturnsS1_low, ReturnsS1, ReturnsS1_high, titS, legend_low, legend_base, legend_high);

clear; clc;


text:
function a = real(a)
%REAL         Real part of gradients
%
%   c = real(a)
%
%

% written  10/16/98     S.M. Rump
% modified 12/18/02     S.M. Rump  improved performance
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  a.x = real(a.x);
  a.dx = real(a.dx);

repo: OpenMAS


text:
function [ p, isgn, ncycle ] = perm_cycle ( n, iopt, p )

%*****************************************************************************80
%
%% PERM_CYCLE analyzes a permutation.
%
%  Discussion:
%
%    The routine will count cycles, find the sign of a permutation,
%    and tag a permutation.
%
%  Example:
%
%    Input:
%
%      N = 9
%      IOPT = 1
%      P = 2, 3, 9, 6, 7, 8, 5, 4, 1
%
%    Output:
%
%      NCYCLE = 3
%      ISGN = +1
%      P_NEW = -2, 3, 9, -6, -7, 8, 5, 4, 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2009
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    A Nijenhuis and H Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the number of objects being permuted.
%
%    Input, integer IOPT, requests tagging.
%    0, the permutation will not be tagged.
%    1, the permutation will be tagged.
%
%    Input, integer P(N), the permutation to be analyzed.
%
%    Output, integer P(N).  If IOPT = 0, then P will be the same as P.
%    If IOPT = 1, then on output, P will be "tagged".
%    That is, one element of every cycle in P will be negated.  In this way,
%    a user can traverse a cycle by starting at any entry I1 of P
%    which is negative, moving to I2 = ABS(P(I1)), then to
%    P(I2), and so on, until returning to I1.
%
%    Output, integer ISGN, the "sign" of the permutation, which is
%    +1 if the permutation is even, -1 if odd.  Every permutation
%    may be produced by a certain number of pairwise switches.
%    If the number of switches is even, the permutation itself is
%    called even.
%
%    Output, integer NCYCLE, the number of cycles in the permutation.
%
  base = 1;

  ierror = perm_check ( n, p, base );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PERM_CYCLE - Fatal error!\n' );
    fprintf ( 1, '  The input array does not represent\n' );
    fprintf ( 1, '  a proper permutation.  In particular, the\n' );
    fprintf ( 1, '  array is missing the value %d\n', ierror );
    error ( 'PERM_CYCLE - Fatal error!' );
  end

  is = 1;
  ncycle = n;

  for i = 1 : n

    i1 = p(i);

    while ( i < i1 )
      ncycle = ncycle - 1;
      i2 = p(i1);
      p(i1) = -i2;
      i1 = i2;
    end

    if ( iopt ~= 0 )
      is = - i4_sign ( p(i) );
    end

    p(i) = i4_sign ( is ) * p(i)

  end

  isgn = 1 - 2 * mod ( n - ncycle, 2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% process text
% propose images
clc;
clear all;
source_file = 'D:\Graduate\Comparison\Comparison_IAFusion\test_Modify\test_name\TNO.xlsx';
source_dir = 'D:\Graduate\Comparison\Comparison_IAFusion\Dataset\New_Test_vi';
save_dir = 'D:\Graduate\Comparison\Comparison_IAFusion\Dataset\PIAFusion\TNO\vi';
if exist(save_dir,'dir')==0.
	 mkdir(save_dir);
end
[NUM,TXT,RAW]=xlsread(source_file);
for i=1:size(NUM, 1)
    if NUM(i) < 10
        image_name = fullfile(source_dir, strcat('0', num2str(NUM(i)), '.bmp'));
        save_name = fullfile(save_dir, strcat('0', num2str(NUM(i)), '.png'));
    elseif NUM(i)>10 & NUM(i) <100
        image_name = fullfile(source_dir, strcat(num2str(NUM(i)), '.bmp'));
        save_name = fullfile(save_dir, strcat(num2str(NUM(i)), '.png'));
    elseif NUM(i)>100
        image_name = fullfile(source_dir, strcat(num2str(NUM(i)), '.bmp'));
        save_name = fullfile(save_dir, strcat(num2str(NUM(i)), '.png'));
    end
%     image_name = fullfile(source_dir, cell2mat(TXT(i)));
%     save_name = fullfile(save_dir, cell2mat(TXT(i)));
    image = imread(image_name);
    [h, w, c] = size(image);
    if c > 1
        image = rgb2gray(image);
    end
    h1 = floor(h / 4) * 4;
    w1 = floor(w / 4) * 4;
    image = image(1:h1, 1:w1);
    imwrite(image, save_name);
    disp(save_name)
end

author: Linfeng-Tang

repo: Image-Fusion


text:
function transform(obj, A, t, varargin)
% TRANSFORM Transform image exterior orientation.

% Input parsing ----------------------------------------------------------------

p = inputParser;
p.addRequired(  'A'              , @(x) isnumeric(x) && size(x,1)==3 && size(x,2)==3);
p.addRequired(  't'              , @(x) isnumeric(x) && numel(x)==3);
% p.addParameter('RedPoi', [0 0 0], @(x) numel(x)==3); % not supported yet!!!
p.parse(A, t, varargin{:});
p = p.Results;
% Clear required inputs to avoid confusion
clear A t

% Start ------------------------------------------------------------------------

% procHierarchy = {'IMG' 'TRANSFORM'};
% msg('S', procHierarchy);
% msg('I', procHierarchy, sprintf('Image label = ''%s''', obj.label));

% Transform --------------------------------------------------------------------

% Unit vectors
xUnitVector = obj.R(:,1);
yUnitVector = obj.R(:,2);
zUnitVector = obj.R(:,3);

% Point cloud for transformation
X = [obj.X0                  obj.Y0                  obj.Z0                    % projection center
     obj.X0 + xUnitVector(1) obj.Y0 + xUnitVector(2) obj.Z0 + xUnitVector(3)   % projection center + xUnitVector
     obj.X0 + yUnitVector(1) obj.Y0 + yUnitVector(2) obj.Z0 + yUnitVector(3)   % projection center + yUnitVector
     obj.X0 + zUnitVector(1) obj.Y0 + zUnitVector(2) obj.Z0 + zUnitVector(3)]; % projection center + zUnitVector

% Homogeneous transformation matrix
H = homotrafo(1, p.A, p.t);

% Homogeneous coordinates
Xh = homocoord(X);

% Transformation of coordinates!
XhTrafo = H * Xh';

% Transformed non homogeneous coordinates
XNew = homocoord(XhTrafo');

% New projection center
X0New = XNew(1,1);
Y0New = XNew(1,2);
Z0New = XNew(1,3);

% New unit vectors
xUnitVectorNew = [XNew(2,1)-X0New XNew(2,2)-Y0New XNew(2,3)-Z0New]';
yUnitVectorNew = [XNew(3,1)-X0New XNew(3,2)-Y0New XNew(3,3)-Z0New]';
zUnitVectorNew = [XNew(4,1)-X0New XNew(4,2)-Y0New XNew(4,3)-Z0New]';

% New rotation matrix
RNew = [xUnitVectorNew yUnitVectorNew zUnitVectorNew];

% New rotation angles
[omeNew, phiNew, kapNew] = R2opk(RNew);

% Save to object
obj.X0  = XNew(1,1);
obj.Y0  = XNew(1,2);
obj.Z0  = XNew(1,3);
obj.ome = omeNew;
obj.phi = phiNew;
obj.kap = kapNew;

% End --------------------------------------------------------------------------

% msg('E', procHierarchy);

end

author: pglira


text:
function stroud_test40 ( )

%*****************************************************************************80
%
%% TEST40 tests TETRA_UNIT_SET, TETRA_SUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_3D_INDEX;

  max_rule = 7;

  x = [ 1.0, 4.0, 1.0, 1.0 ];
  y = [ 2.0, 2.0, 3.0, 2.0 ];
  z = [ 6.0, 6.0, 6.0, 8.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST40\n' );
  fprintf ( 1, '  TETRA_UNIT_SET sets quadrature rules\n' );
  fprintf ( 1, '    for the unit tetrahedron;\n' );
  fprintf ( 1, '  TETRA_SUM applies them to an arbitrary tetrahedron.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tetrahedron vertices:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 4
    fprintf ( 1, '  %3f  %3f  %3f\n', x(i), y(i), z(i) );
  end

  for ilo = 1 : 5 : max_rule

    ihi = min ( ilo +  4, max_rule );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Rule:    ' );
    for rule = ilo : ihi
      fprintf ( 1, '%7d       ', rule );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Function\n' );
    fprintf ( 1, '\n' );

    num = function_3d_num ( );

    for i = 1 : num

      FUNC_3D_INDEX = i;

      for rule = ilo : ihi

        order = tetra_unit_size ( rule );

        [ xtab, ytab, ztab, weight ] = tetra_unit_set ( rule, order );

        result(rule) = tetra_sum ( 'function_3d', x, y, z, order, xtab, ytab, ztab, weight );

      end

      fname = function_3d_name ( i );

      fprintf ( 1, '  %s', fname );
      for rule = ilo : ihi
        fprintf ( 1, '  %12f', result(rule) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function plotData(x, y)
%PLOTDATA Plots the data points x and y into a new figure
%   PLOTDATA(x,y) plots the data points and gives the figure axes labels of
%   population and profit.

figure; % open a new figure window
plot(x, y, 'rx', 'MarkerSize', 10);
ylabel('Profit in $10,000s');
xlabel('Population of City in 10,000s');

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the training data into a figure using the
%               "figure" and "plot" commands. Set the axes labels using
%               the "xlabel" and "ylabel" commands. Assume the
%               population and revenue data have been passed in
%               as the x and y arguments of this function.
%
% Hint: You can use the 'rx' option with plot to have the markers
%       appear as red crosses. Furthermore, you can make the
%       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

% ============================================================

end

author: loserChen

repo: Coursera-MachineLearning


text:
function [newModel, rxnToCheck]=expandModel(model)
% expandModel
%   Expands a model which uses several gene associations for one reaction.
%   Each such reaction is split into several reactions, each under the control
%   of only one gene.
%
% Input:
%   model       model structure
%
% Output:
%   newModel    model structure with separate reactions for iso-enzymes, where
%               the reaction ids are renamed as to id_EXP_1, id_EXP_2, etc.
%   rxnToCheck  cell array with original reaction identifiers for those
%               that contained nested and/or-relationships in grRules.
%
%   NOTE: grRules strings that involve nested expressions of 'and' and 'or'
%         might not be parsed correctly if they are not standardized (if the
%         standardizeGrRules functions was not first run on the model). For
%         those reactions, it is therefore advisable to inspect the reactions in
%         rxnToCheck to confirm correct model expansion.
%
%   Usage: [newModel, rxnToCheck]=expandModel(model)

%Check how many reactions we will create (the number of or:s in the GPRs).
%This way, we can preallocate all fields and save much computation time

numOrs = count(model.grRules, ' or ');
toAdd = sum(numOrs);
prevNumRxns = length(model.rxns);
rxnToCheck={};
if toAdd > 0
    %Calculate indices to expand
    %For example, if a reaction with index x has 2 or:s, meaning it has 3
    %reactions after the split, we should add two copies of this reaction
    %For fields that should just be copied to the new reactions, we just keep
    %track of that there are two copies, i.e., we add x x to this vector.
    %That is exactly what repelem does for us.
    cpyIndices = repelem(1:prevNumRxns, numOrs);

    %Copy all fields that should just be copied
    model.S=[model.S model.S(:,cpyIndices)];
    model.rxnNames=[model.rxnNames;model.rxnNames(cpyIndices)];
    model.lb=[model.lb;model.lb(cpyIndices)];
    model.ub=[model.ub;model.ub(cpyIndices)];
    model.rev=[model.rev;model.rev(cpyIndices)];
    model.c=[model.c;model.c(cpyIndices)];
    if isfield(model,'subSystems')
        model.subSystems=[model.subSystems;model.subSystems(cpyIndices)];
    end
    if isfield(model,'eccodes')
        model.eccodes=[model.eccodes;model.eccodes(cpyIndices)];
    end
    if isfield(model,'equations')
        model.equations=[model.equations;model.equations(cpyIndices)];
    end
    if isfield(model,'rxnMiriams')
        model.rxnMiriams=[model.rxnMiriams;model.rxnMiriams(cpyIndices)];
    end
    if isfield(model,'rxnComps')
        model.rxnComps=[model.rxnComps;model.rxnComps(cpyIndices)];
    end
    if isfield(model,'rxnFrom')
        model.rxnFrom=[model.rxnFrom;model.rxnFrom(cpyIndices)];
    end
    if isfield(model,'rxnNotes')
        model.rxnNotes=[model.rxnNotes;model.rxnNotes(cpyIndices)];
    end
    if isfield(model,'rxnReferences')
        model.rxnReferences=[model.rxnReferences;model.rxnReferences(cpyIndices)];
    end
    if isfield(model,'rxnConfidenceScores')
        model.rxnConfidenceScores=[model.rxnConfidenceScores;model.rxnConfidenceScores(cpyIndices)];
    end

    %now expand the more complex fields - will be filled in later
    model.rxns=[model.rxns;cell(toAdd,1)];
    model.grRules=[model.grRules;cell(toAdd,1)];
    model.rxnGeneMat=[model.rxnGeneMat;sparse(toAdd,size(model.rxnGeneMat,2))];

    %Loop throught those reactions and fill in the expanded data
    nextIndex = prevNumRxns + 1;
    for i=1:prevNumRxns
        if (numOrs(i) > 0)
            %Check that it doesn't contain nested 'and' and 'or' relations and
            %print a warning if it does
            if ~isempty(strfind(model.grRules{i},' and '))
                rxnToCheck{end+1,1}=model.rxns{i};
            end

            %Get rid of all '(' and ')' since I'm not looking at complex stuff
            %anyways
            geneString=model.grRules{i};
            geneString=strrep(geneString,'(','');
            geneString=strrep(geneString,')','');
            geneString=strrep(geneString,' or ',';');

            %Split the string into gene names
            geneNames=regexp(geneString,';','split');

            %Update the reaction to only use the first gene
            model.grRules{i}=['(' geneNames{1} ')'];
            %Find the gene in the gene list If ' and ' relationship, first
            %split the genes
            model.rxnGeneMat(i,:)=0;
            if ~isempty(strfind(geneNames(1),' and '))
                andGenes=regexp(geneNames{1},' and ','split');
                model.rxnGeneMat(i,ismember(model.genes,andGenes)) = 1;
            else
                [~, index]=ismember(geneNames(1),model.genes);
                model.rxnGeneMat(i,index)=1;
            end

            %Insert the reactions at the end of the model and without
            %allocating space. This is not nice, but ok for now
            for j=2:numel(geneNames)
                ind = nextIndex+j-2;
                model.rxns{ind}=[model.rxns{i} '_EXP_' num2str(j)];

                model.grRules{ind}=['(' geneNames{j} ')'];

                if ~isempty(strfind(geneNames(j),' and '))
                    andGenes=regexp(geneNames{j},' and ','split');
                    model.rxnGeneMat(ind,ismember(model.genes,andGenes)) = 1;
                else
                    model.rxnGeneMat(ind,ismember(model.genes,geneNames(j))) = 1;
                end
            end
            model.rxns{i}=[model.rxns{i}, '_EXP_1'];
            nextIndex = nextIndex + numOrs(i);
        end
    end
    newModel=model;
else
    %There are no reactions to expand, return the model as is
    newModel=model;
end

%Fix grRules and reconstruct rxnGeneMat
[grRules,rxnGeneMat] = standardizeGrRules(newModel,true);
newModel.grRules = grRules;
newModel.rxnGeneMat = rxnGeneMat;
end

author: SysBioChalmers

repo: RAVEN


text:
function XCorr=SetNaNs1(X,ExMin,ExMax,EmMin,EmMax,WidthMin,WidthMax,Value,MaxOrder,Mode);
%SETNANS1 fluorescence artifact treatment
%
%XCorr=SetNaNs1(X,ExMin,ExMax,EmMin,EmMax,WidthMin,WidthMax,Value,MaxOrder,Mode);
%
%This algorithm takes a fluorescence landscape with
%excitation ranges and emission ranges defined by
%ExMin, ExMax, EmMin and EmMax. Given these specifications
%the algorithm identifies the Rayleigh-peak and replaces
%a certain number of neighboring wavelengths at each side
%of the peak with the 'Value' multiplied by
%the existing value in the range, - use NaN to set missing.
%
%Generally the peaks broaden as the excitation wavelength increase.
%The 'WidthMin' is the width of the Rayleigh peaks at the lowest
%excitation wavelength and 'WidthMax' is the peakwidth at the highest
%excitation wavelength. Giv also the width of the 2nd order Rayleigh
%peak of the first spectrum as 'Width2nd' since 2nd order scattering may be
%much wider than first order. The algorithm replaces up to and including
%scattering of 'MaxOrder' order starting from the 1st (no 0th order).
%'Mode'=1 uses the lowest emission wavelength as low limit, hence
%removing all emission wavelengths on the low side of the 0th order
%Rayleigh peak. 'Mode'=2 removes only the peak itself.
%
%The landscape is contained in the matrix X and each row in the matrix
%represents one emission spectrum. Hence, X has dimensions (a x b) where
%a is the number of excitation lines, and b is the number
%of measured emission wavelengths.
%
%                                  b
%          -------------------------
%          | 1st emission spectrum |
%          | 2nd emission spectrum |
%          |         ....          |
%          | ath emission spectrum |
%        a -------------------------
%
% %XCorr=SetNaNs1(X,ExMin,ExMax,EmMin,EmMax,WidthMin,WidthMax,Value,MaxOrder,Mode);

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

XCorr=X;
[a b]=size(X);

exax=ExMin:(ExMax-ExMin)/(a-1):ExMax;
emax=EmMin:(EmMax-EmMin)/(b-1):EmMax;

if WidthMin==WidthMax,
   Width=WidthMin*ones(1,a);
else
   Width=WidthMin:(WidthMax-WidthMin)/(a-1):WidthMax;
end;

WFactor=[1 1.4 1.7 1.9 2];

for order=1:MaxOrder,

   for excount=1:a,

      exc=order*exax(excount);

      W = WFactor(order)*Width(excount);

      rangemin = max(find(emax<=(exc-W)));
      rangemax = min(find(emax>=(exc+W)));

      if (exc-W) < EmMin,
         rangemin = 1;
      end;

      if (exc-W) > EmMax,
         rangemin = 0;
      end;

      if (exc+W) < EmMin,
         rangemax = 0;
      end;

      if (exc+W) > EmMax,
         rangemax = b;
      end;

      if ~(rangemax*rangemin==0),
         if Mode==1 & order==1;
            rangemin = 1;
         end;
         XCorr(excount,rangemin:rangemax)=Value*X(excount,rangemin:rangemax);
         [rangemin rangemax];
      end;

   end;
end;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function volume = polyhedron_volume_3d_2 ( coord, maxorder, face_num, node, ...
  node_num, order )

%*****************************************************************************80
%
%% POLYHEDRON_VOLUME_3D_2 computes the volume of a polyhedron in 3D.
%
%  Discussion:
%
%    The computation is not valid unless the faces of the polyhedron
%    are planar polygons.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Allen Van Gelder,
%    Efficient Computation of Polygon Area and Polyhedron Volume,
%    Graphics Gems V, edited by Alan Paeth,
%    AP Professional, 1995.
%
%  Parameters:
%
%    Input, real COORD(3,NODE_NUM), the vertices.
%    The vertices may be listed in any order.
%
%    Input, integer MAXORDER, the maximum number of vertices that make
%    up a face of the polyhedron.
%
%    Input, integer FACE_NUM, the number of faces of the polyhedron.
%
%    Input, integer NODE(FACE_NUM,MAXORDER).  Face I is defined by
%    the vertices NODE(I,1) through NODE(I,ORDER(I)).  These vertices
%    are listed in neighboring order.
%
%    Input, integer NODE_NUM, the number of points stored in COORD.
%
%    Input, integer ORDER(FACE_NUM), the number of vertices making up
%    each face.
%
%    Output, real VOLUME, the volume of the polyhedron.
%
  dim_num = 3;

  volume = 0.0;

  for face = 1 : face_num

    v(1:dim_num) = 0.0;
%
%  Compute the area vector for this face.
%
    for j = 1 : order(face)

      k1 = node(face,j);

      if ( j < order(face) )
        k2 = node(face,j+1);
      else
        k2 = node(face,1);
      end
%
%  Compute the cross product.
%
      normal(1) = coord(2,k1) * coord(3,k2) - coord(3,k1) * coord(2,k2);
      normal(2) = coord(3,k1) * coord(1,k2) - coord(1,k1) * coord(3,k2);
      normal(3) = coord(1,k1) * coord(2,k2) - coord(2,k1) * coord(1,k2);

      v(1:dim_num) = v(1:dim_num) + normal(1:dim_num);

    end
%
%  Area vector dot any vertex.
%
    k = node(face,1);
    volume = volume + v(1:dim_num) * coord(1:dim_num,k);

  end

  volume = volume / 6.0;

  return
end


text:
function SO3VF = rdivide(SO3VF1, SO3VF2)
% overloads |SO3VF1 ./ SO3VF2|
%
% Syntax
%   SO3VF = SO3VF1 ./ SO3VF2
%   SO3VF = SO3VF1 ./ a
%   SO3VF = a ./ SO3VF2
%   SO3VF = SO3F .* SO3VF1
%   SO3VF = SO3VF1 .* SO3F
%
% Input
%  SO3VF1, SO3VF2 - @SO3VectorField
%  a - double, @vector3d
%  SO3F - @SO3Fun
%
% Output
%  SO3VF - @SO3VectorField
%

if isnumeric(SO3VF1)
  SO3VF = SO3VectorFieldHandle(@(rot) SO3VF1./ SO3VF2.eval(rot),SO3VF2.SRight,SO3VF2.SLeft);
  return
end
if isa(SO3VF2,'SO3Fun')   % prevent quadrature in case of SO3FunHarmonic
  SO3VF = SO3VectorFieldHandle(@(rot) SO3VF1.eval(rot)./ SO3VF2.eval(rot),SO3VF2.SRight,SO3VF2.SLeft);
  return
end
SO3VF = times(SO3VF1,1./SO3VF2);

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
% Run sanity checks on a given set of image pairs.
function [sanityCheckPassed, sanityCheckFailed, failedChecks] = performSanityCheck(imagePairs, maxSanityCheckFails)

    numWorkers = length(imagePairs);
    sanityCheckPassed = [];
    sanityCheckFailed = [];
    for workerIdx = 1:numWorkers

        [sanityCheck, failedChecks(:,workerIdx)] = performSanityCheckForWorker(imagePairs{workerIdx}, maxSanityCheckFails);
        if sanityCheck == true
            sanityCheckPassed = cat(1, sanityCheckPassed, workerIdx);
        else
            sanityCheckFailed = cat(1, sanityCheckFailed, workerIdx);
        end

    end

function [sanityCheckPassed, numChecksFailed] = performSanityCheckForWorker(imagePairs, maxSanityCheckFails)

    numChecksFailed = zeros(length(imagePairs), 1);
    sanityCheckPassed = true;
    for imagePairIdx = 1:length(imagePairs)

        % Check if the current pair is used as a sanity check.
        if imagePairs(imagePairIdx).isSanityCheck == true

            % Check if the current pair passes the sanity check.
            if ~isempty(strfind(imagePairs(imagePairIdx).firstImage.name, '_good')) && imagePairs(imagePairIdx).vote == 2
                % It failed the sanity check.
                numChecksFailed(imagePairIdx) = 1;
            end
            if ~isempty(strfind(imagePairs(imagePairIdx).firstImage.name, '_bad')) && imagePairs(imagePairIdx).vote == 1
                % It failed the sanity check.
                numChecksFailed(imagePairIdx) = 1;
            end

        end

    end

    if sum(numChecksFailed) > maxSanityCheckFails
        sanityCheckPassed = false;
    end

author: thomas-koehler

repo: SupER


text:
% Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, see <http://www.gnu.org/licenses/>.
%
% usage: w = gausswin(L, a)
%
% Generate an L-point gaussian window of the given width. Use larger a
% for a narrow window.  Use larger L for a smoother curve.
%
%     w = exp ( -(a*x)^2/2 )
%
% for x = linspace(-(L-1)/L, (L-1)/L, L)

function x = gausswin(L, w)

  if nargin < 1 || nargin > 2
    help(mfilename);
  end
  if nargin == 1, w = 2.5; end % if

  if L==1
    x = 1;
  else
    x = exp ( -0.5 * ( w/(L-1) * ( -(L-1) : 2 : L-1 )' ) .^ 2 );
  end

end


text:
function git(varargin)
%   A thin MATLAB wrapper for the git source control system
%   The short instructions are:
%       Use this exactly as you would use the OS command-line Git.
%
%   The long instructions are:
%       This is not meant to be a comprehensive guide to near-omnipotent
%       Git source control. That is here:
%           http://git-scm.com/documentation
%
%       Common MATLAB workflow:
%
%       % Creates initial repository tracking all files under some root
%       % folder
%       >> cd ~/
%       >> git init
%
%       % Shows changes made to all files in repo (none so far)
%       >> git status
%
%       % Create a new file and add some code
%       >> edit foo.m
%
%       % Check repo status, after new file created
%       >> git status
%
%       % Add foo.m to your repo
%       >> git add foo.m
%
%       % Commit your changes to a new branch, with comments
%       >> git commit -m 'Created new file, foo.m'
%
%       % Other useful commands (replace ellipses with appropriate args)
%       >> git checkout ...       % To return to an earlier node in branch
%       >> git branch ...         % To create or move to another branch
%       >> git merge ...          % When merge conflicts arise
%       >> git diff ...           % See line-by-line changes
%
%   Useful resources:
%       1. GitX: A visual interface for Git on the OS X client
%       2. Github.com: Remote hosting for Git repos
%       3. Git on Wikipedia: Further reading
%
% v0.1,     27 October, 2010 -- Initial support for OS X & Linux, untested
%                               on PCs, but expected to work
%
% Author: Manu Raghavan
% Free for commercial use

% if ismac || isunix
%     [status] = system('git');
%     assert(status==1, ...
%         strcat('git is not installed\n',...
%                'Download it at http://git-scm.com/download'));
% elseif ispc
%     warning('Make sure git is installed');
% end

arguments = parse(varargin{:});
[smaz,result] = system(['git ',arguments]);
disp(result)
end

function space_delimited_list = parse(varargin)
    space_delimited_list = cell2mat(...
                cellfun(@(s)([s,' ']),varargin,'UniformOutput',false));
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [mh] = spm_mci_popdef (scale,tune,samp)
% Set default parameters for population MCMC
% FORMAT [mh] = spm_mci_popdef (scale,tune,samp)
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_popdef.m 6548 2015-09-11 12:39:47Z will $

mh.nscale=scale;
mh.ntune=tune;
mh.nsamp=samp;
mh.ind_samp=[scale+tune+1:scale+tune+samp];
mh.J=1; % Number of temperatures
mh.gprob=0;
mh.remove_burn_in=0;
mh.verbose=0;

author: spm


text:
function varargout = log10(varargin)

varargout{1} = log(varargin{1})/log(10);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function OUT = Date2Num(dates,frequency)
% =========================================================================
% Convert a cell array of dates into a double array. The first period of
% the year is coded as a round number. For quarterly data the step change
% is 0.25, so that 1984Q2 => 1984.25; For  monthly data the step change
% is 0.833, so that 1984M2 => 1984.0833.
% =======================================================================
% OUT = Date2Num(DATA)
% -----------------------------------------------------------------------
% INPUTS
%	- DATA: a (T x 1) vector with dates in cell format
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - frequency : monthly ('m'), quarterly ('q') [dflt] or annual ('y')
% -----------------------------------------------------------------------
% OUTPUT
%	- OUT:  a (T x 1) vector with dates in numeric format
% -----------------------------------------------------------------------
% EXAMPLE
% 	dates = {'1992M1';'1992M2';'1992M3';'1992M4';'1992M5';};
%   OUT = Date2Num(dates,'m')
% -----------------------------------------------------------------------
% RELATED
%   - DatesCreate, DatesCount
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2015. Updated November 2020
% -----------------------------------------------------------------------

format long
% Check inputs
if ~exist('frequency','var')
    frequency = 'q';
end
[n,c] = size(dates);

% Convert cell array into double array
if strcmp(frequency,'q')
    OUT = nan(n,c);
    for i=1:n
        for j=1:c
            aux = dates{i,j};
            year = str2double(aux(1:4));
            period = str2double(aux(6));
            OUT(i,j) = year + period/4 -0.25;
        end
    end
elseif strcmp(frequency,'m')
    OUT = nan(n,c);
    for i=1:n
        for j=1:c
            aux = dates{i,j};
            year = str2double(aux(1:4));
            period = str2double(aux(6:end));
            OUT(i,j) = year + period/12 -0.08333333333333333333333333333333;
        end
    end
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function VO = spm_write_sn(V,prm,flags,extras)
% Write Out Warped Images.
% FORMAT VO = spm_write_sn(V,matname,flags,msk)
% V         - Images to transform (filenames or volume structure).
% matname   - Transformation information (filename or structure).
% flags     - flags structure, with fields...
%           interp   - interpolation method (0-7)
%           wrap     - wrap edges (e.g., [1 1 0] for 2D MRI sequences)
%           vox      - voxel sizes (3 element vector - in mm)
%                      Non-finite values mean use template vox.
%           bb       - bounding box (2x3 matrix - in mm)
%                      Non-finite values mean use template bb.
%           preserve - either 0 or 1.  A value of 1 will "modulate"
%                      the spatially normalised images so that total
%                      units are preserved, rather than just
%                      concentrations.
%           prefix   - Prefix for normalised images. Defaults to 'w'.
% msk       - An optional cell array for masking the spatially
%             normalised images (see below).
%
% Warped images are written prefixed by "w".
%
% Non-finite vox or bounding box suggests that values should be derived
% from the template image.
%
% Don't use interpolation methods greater than one for data containing
% NaNs.
% _______________________________________________________________________
%
% FORMAT msk = spm_write_sn(V,matname,flags,'mask')
% V         - Images to transform (filenames or volume structure).
% matname   - Transformation information (filename or structure).
% flags     - flags structure, with fields...
%           wrap     - wrap edges (e.g., [1 1 0] for 2D MRI sequences)
%           vox      - voxel sizes (3 element vector - in mm)
%                      Non-finite values mean use template vox.
%           bb       - bounding box (2x3 matrix - in mm)
%                      Non-finite values mean use template bb.
% msk       - a cell array for masking a series of spatially normalised
%             images.
%
%
% _______________________________________________________________________
%
% FORMAT VO = spm_write_sn(V,prm,'modulate')
% V         - Spatially normalised images to modulate (filenames or
%             volume structure).
% prm       - Transformation information (filename or structure).
%
%  After nonlinear spatial normalization, the relative volumes of some
%  brain structures will have decreased, whereas others will increase.
%  The resampling of the images preserves the concentration of pixel
%  units in the images, so the total counts from structures that have
%  reduced volumes after spatial normalization will be reduced by an
%  amount proportional to the volume reduction.
%
%  This routine rescales images after spatial normalization, so that
%  the total counts from any structure are preserved.  It was written
%  as an optional step in performing voxel based morphometry.
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

if isempty(V), return; end;

if ischar(prm), prm = load(prm);  end;
if ischar(V),   V   = spm_vol(V); end;

if nargin==3 && ischar(flags) && strcmpi(flags,'modulate'),
    if nargout==0,
        modulate(V,prm);
    else
        VO = modulate(V,prm);
    end;
    return;
end;

def_flags = struct('interp',1,'vox',NaN,'bb',NaN,'wrap',[0 0 0],'preserve',0,...
                   'prefix','w');

if nargin < 3,
    flags = def_flags;
else
    fnms = fieldnames(def_flags);
    for i=1:length(fnms),
        if ~isfield(flags,fnms{i}),
            flags.(fnms{i}) = def_flags.(fnms{i});
        end;
    end;
end;

[x,y,z,mat] = get_xyzmat(prm,flags.bb,flags.vox);

if nargin==4,
    if ischar(extras) && strcmpi(extras,'mask'),
        VO = get_snmask(V,prm,x,y,z,flags.wrap);
        return;
    end;
    if iscell(extras),
        msk = extras;
    end;
end;

if nargout>0 && length(V)>8,
    error('Too many images to save in memory');
end;

if ~exist('msk','var')
    msk         = get_snmask(V,prm,x,y,z,flags.wrap);
end;

if nargout==0,
    if isempty(prm.Tr),
        affine_transform(V,prm,x,y,z,mat,flags,msk);
    else
        nonlin_transform(V,prm,x,y,z,mat,flags,msk);
    end;
else
    if isempty(prm.Tr),
        VO = affine_transform(V,prm,x,y,z,mat,flags,msk);
    else
        VO = nonlin_transform(V,prm,x,y,z,mat,flags,msk);
    end;
end;

return;
%_______________________________________________________________________

%_______________________________________________________________________
function VO = affine_transform(V,prm,x,y,z,mat,flags,msk)

[X,Y] = ndgrid(x,y);
d     = [flags.interp*[1 1 1]' flags.wrap(:)];

spm_progress_bar('Init',numel(V),'Resampling','volumes/slices completed');
for i=1:numel(V),
    VO     = make_hdr_struct(V(i),x,y,z,mat, flags.prefix);
    if flags.preserve
        VO.fname = prepend(VO.fname,'m');
    end
    detAff = det(prm.VF(1).mat*prm.Affine/prm.VG(1).mat);
    if flags.preserve, VO.pinfo(1:2,:) = VO.pinfo(1:2,:)/detAff; end;

    %Dat= zeros(VO.dim(1:3));
    Dat = single(0);
    Dat(VO.dim(1),VO.dim(2),VO.dim(3)) = 0;

    C = spm_bsplinc(V(i),d);

    for j=1:length(z),   % Cycle over planes
        [X2,Y2,Z2]  = mmult(X,Y,z(j),V(i).mat\prm.VF(1).mat*prm.Affine);
        dat         = spm_bsplins(C,X2,Y2,Z2,d);
        if flags.preserve, dat = dat*detAff; end;
        dat(msk{j}) = NaN;

        Dat(:,:,j) = single(dat);

        if numel(V)<5, spm_progress_bar('Set',i-1+j/length(z)); end;
    end;
    if nargout~=0,
        VO.pinfo  = [1 0]';
        VO.dt     = [spm_type('float32') spm_platform('bigend')];
        VO.dat    = Dat;
    else
        spm_write_vol(VO, Dat);
    end;
    spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');
return;
%_______________________________________________________________________

%_______________________________________________________________________
function VO = nonlin_transform(V,prm,x,y,z,mat,flags,msk)

[X,Y] = ndgrid(x,y);
Tr = prm.Tr;
BX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1);
BY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1);
BZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1);
if flags.preserve,
    DX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1,'diff');
    DY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1,'diff');
    DZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1,'diff');
end;
d  = [flags.interp*[1 1 1]' flags.wrap(:)];

spm_progress_bar('Init',numel(V),'Resampling','volumes completed');
for i=1:numel(V),
    VO     = make_hdr_struct(V(i),x,y,z,mat, flags.prefix);
    if flags.preserve
        VO.fname = prepend(VO.fname,'m');
    end
    detAff = det(prm.VF(1).mat*prm.Affine/prm.VG(1).mat);

    % Accumulate data
    %Dat= zeros(VO.dim(1:3));
    Dat = single(0);
    Dat(VO.dim(1),VO.dim(2),VO.dim(3)) = 0;

    C = spm_bsplinc(V(i),d);

    for j=1:length(z),   % Cycle over planes
        % Nonlinear deformations
        %----------------------------------------------------------------------------
        tx = get_2Dtrans(Tr(:,:,:,1),BZ,j);
        ty = get_2Dtrans(Tr(:,:,:,2),BZ,j);
        tz = get_2Dtrans(Tr(:,:,:,3),BZ,j);
        X1 = X    + BX*tx*BY';
        Y1 = Y    + BX*ty*BY';
        Z1 = z(j) + BX*tz*BY';

        [X2,Y2,Z2]  = mmult(X1,Y1,Z1,V(i).mat\prm.VF(1).mat*prm.Affine);
        dat         = spm_bsplins(C,X2,Y2,Z2,d);
        dat(msk{j}) = NaN;

        if ~flags.preserve,
            Dat(:,:,j) = single(dat);
        else
            j11 = DX*tx*BY' + 1; j12 = BX*tx*DY';     j13 = BX*get_2Dtrans(Tr(:,:,:,1),DZ,j)*BY';
            j21 = DX*ty*BY';     j22 = BX*ty*DY' + 1; j23 = BX*get_2Dtrans(Tr(:,:,:,2),DZ,j)*BY';
            j31 = DX*tz*BY';     j32 = BX*tz*DY';     j33 = BX*get_2Dtrans(Tr(:,:,:,3),DZ,j)*BY' + 1;

            % The determinant of the Jacobian reflects relative volume changes.
            %------------------------------------------------------------------
            dat       = dat .* (j11.*(j22.*j33-j23.*j32) - j21.*(j12.*j33-j13.*j32) + j31.*(j12.*j23-j13.*j22)) * detAff;
            Dat(:,:,j) = single(dat);
        end;
        if numel(V)<5, spm_progress_bar('Set',i-1+j/length(z)); end;
    end;
    if nargout==0,
        if flags.preserve, VO = rmfield(VO,'pinfo'); end
        VO = spm_write_vol(VO,Dat);
    else
        VO.pinfo  = [1 0]';
        VO.dt     = [spm_type('float32') spm_platform('bigend')];
        VO.dat    = Dat;
    end;
    spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');
return;
%_______________________________________________________________________

%_______________________________________________________________________
function VO = modulate(V,prm)

spm_progress_bar('Init',numel(V),'Modulating','volumes completed');
for i=1:numel(V),
    VO          = V(i);
    VO          = rmfield(VO,'pinfo');
    VO.fname    = prepend(VO.fname,'m');
    detAff      = det(prm.VF(1).mat*prm.Affine/prm.VG(1).mat);
    %Dat        = zeros(VO.dim(1:3));
    Dat         = single(0);
    Dat(VO.dim(1),VO.dim(2),VO.dim(3)) = 0;
    [x,y,z,mat] = get_xyzmat(prm,NaN,NaN,VO);

    if sum((mat(:)-VO.mat(:)).^2)>1e-7, error('Orientations not compatible'); end;

    Tr = prm.Tr;

    if isempty(Tr),
        for j=1:length(z),   % Cycle over planes
            dat        = spm_slice_vol(V(i),spm_matrix([0 0 j]),V(i).dim(1:2),0);
            Dat(:,:,j) = single(dat);
            if numel(V)<5, spm_progress_bar('Set',i-1+j/length(z)); end;
        end;
    else
        BX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1);
        BY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1);
        BZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1);
        DX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1,'diff');
        DY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1,'diff');
        DZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1,'diff');

        for j=1:length(z),   % Cycle over planes

            tx = get_2Dtrans(Tr(:,:,:,1),BZ,j);
            ty = get_2Dtrans(Tr(:,:,:,2),BZ,j);
            tz = get_2Dtrans(Tr(:,:,:,3),BZ,j);

            j11 = DX*tx*BY' + 1; j12 = BX*tx*DY';     j13 = BX*get_2Dtrans(Tr(:,:,:,1),DZ,j)*BY';
            j21 = DX*ty*BY';     j22 = BX*ty*DY' + 1; j23 = BX*get_2Dtrans(Tr(:,:,:,2),DZ,j)*BY';
            j31 = DX*tz*BY';     j32 = BX*tz*DY';     j33 = BX*get_2Dtrans(Tr(:,:,:,3),DZ,j)*BY' + 1;

            % The determinant of the Jacobian reflects relative volume changes.
            %------------------------------------------------------------------
            dat        = spm_slice_vol(V(i),spm_matrix([0 0 j]),V(i).dim(1:2),0);
            dat        = dat .* (j11.*(j22.*j33-j23.*j32) - j21.*(j12.*j33-j13.*j32) + j31.*(j12.*j23-j13.*j22)) * detAff;
            Dat(:,:,j) = single(dat);

            if numel(V)<5, spm_progress_bar('Set',i-1+j/length(z)); end;
        end;
    end;

    if nargout==0,
        VO = spm_write_vol(VO,Dat);
    else
        VO.pinfo  = [1 0]';
        VO.dt     = [spm_type('float32') spm_platform('bigend')];
        VO.dat    = Dat;
    end;
    spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');
return;
%_______________________________________________________________________

%_______________________________________________________________________
function VO   = make_hdr_struct(V,x,y,z,mat,prefix)
VO            = V;
VO.fname      = prepend(V.fname,prefix);
VO.mat        = mat;
VO.dim(1:3)   = [length(x) length(y) length(z)];
VO.pinfo      = V.pinfo;
VO.descrip    = 'spm - 3D normalized';
return;
%_______________________________________________________________________

%_______________________________________________________________________
function T2 = get_2Dtrans(T3,B,j)
d   = [size(T3) 1 1 1];
tmp = reshape(T3,d(1)*d(2),d(3));
T2  = reshape(tmp*B(j,:)',d(1),d(2));
return;
%_______________________________________________________________________

%_______________________________________________________________________
function PO = prepend(PI,pre)
[pth,nm,xt,vr] = spm_fileparts(deblank(PI));
PO             = fullfile(pth,[pre nm xt vr]);
return;
%_______________________________________________________________________

%_______________________________________________________________________
function Mask = getmask(X,Y,Z,dim,wrp)
% Find range of slice
tiny = 5e-2;
Mask = true(size(X));
if ~wrp(1), Mask = Mask & (X >= (1-tiny) & X <= (dim(1)+tiny)); end;
if ~wrp(2), Mask = Mask & (Y >= (1-tiny) & Y <= (dim(2)+tiny)); end;
if ~wrp(3), Mask = Mask & (Z >= (1-tiny) & Z <= (dim(3)+tiny)); end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [X2,Y2,Z2] = mmult(X1,Y1,Z1,Mult)
if length(Z1) == 1,
    X2= Mult(1,1)*X1 + Mult(1,2)*Y1 + (Mult(1,3)*Z1 + Mult(1,4));
    Y2= Mult(2,1)*X1 + Mult(2,2)*Y1 + (Mult(2,3)*Z1 + Mult(2,4));
    Z2= Mult(3,1)*X1 + Mult(3,2)*Y1 + (Mult(3,3)*Z1 + Mult(3,4));
else
    X2= Mult(1,1)*X1 + Mult(1,2)*Y1 + Mult(1,3)*Z1 + Mult(1,4);
    Y2= Mult(2,1)*X1 + Mult(2,2)*Y1 + Mult(2,3)*Z1 + Mult(2,4);
    Z2= Mult(3,1)*X1 + Mult(3,2)*Y1 + Mult(3,3)*Z1 + Mult(3,4);
end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function msk = get_snmask(V,prm,x,y,z,wrap)
% Generate a mask for where there is data for all images
%-----------------------------------------------------------------------
msk = cell(length(z),1);
t1 = cat(3,V.mat);
t2 = cat(1,V.dim);
t  = [reshape(t1,[16 length(V)])' t2(:,1:3)];
Tr = prm.Tr;
[X,Y] = ndgrid(x,y);

BX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1);
BY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1);
BZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1);

if numel(V)>1 && any(any(diff(t,1,1))),
    spm_progress_bar('Init',length(z),'Computing available voxels','planes completed');
    for j=1:length(z),   % Cycle over planes
        Count = zeros(length(x),length(y));
        if isempty(Tr),
            % Generate a mask for where there is data for all images
            %----------------------------------------------------------------------------
            for i=1:numel(V),
                [X2,Y2,Z2] = mmult(X,Y,z(j),V(i).mat\prm.VF(1).mat*prm.Affine);
                Count      = Count + getmask(X2,Y2,Z2,V(i).dim(1:3),wrap);
            end;
        else
            % Nonlinear deformations
            %----------------------------------------------------------------------------
            X1 = X    + BX*get_2Dtrans(Tr(:,:,:,1),BZ,j)*BY';
            Y1 = Y    + BX*get_2Dtrans(Tr(:,:,:,2),BZ,j)*BY';
            Z1 = z(j) + BX*get_2Dtrans(Tr(:,:,:,3),BZ,j)*BY';

            % Generate a mask for where there is data for all images
            %----------------------------------------------------------------------------
            for i=1:numel(V),
                [X2,Y2,Z2] = mmult(X1,Y1,Z1,V(i).mat\prm.VF(1).mat*prm.Affine);
                Count      = Count + getmask(X2,Y2,Z2,V(i).dim(1:3),wrap);
            end;
        end;
        msk{j} = uint32(find(Count ~= numel(V)));
        spm_progress_bar('Set',j);
    end;
    spm_progress_bar('Clear');
else
    for j=1:length(z), msk{j} = uint32([]); end;
end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [x,y,z,mat] = get_xyzmat(prm,bb,vox,VG)
% The old voxel size and origin notation is used here.
% This requires that the position and orientation
% of the template is transverse.  It would not be
% straitforward to account for templates that are
% in different orientations because the basis functions
% would no longer be seperable.  The seperable basis
% functions mean that computing the deformation field
% from the parameters is much faster.

% bb  = sort(bb);
% vox = abs(vox);

if nargin<4,
    VG = prm.VG(1);

    if all(~isfinite(bb(:))) && all(~isfinite(vox(:))),
        x   = 1:VG.dim(1);
        y   = 1:VG.dim(2);
        z   = 1:VG.dim(3);
        mat = VG.mat;
        return;
    end
end

[bb0,vox0] = bbvox_from_V(VG);
if ~all(isfinite(vox(:))), vox = vox0; end;
if ~all(isfinite(bb(:))),  bb  = bb0;  end;

msk       = find(vox<0);
bb        = sort(bb);
bb(:,msk) = flipud(bb(:,msk));

% Adjust bounding box slightly - so it rounds to closest voxel.
% Comment out if not needed.
%bb(:,1) = round(bb(:,1)/vox(1))*vox(1);
%bb(:,2) = round(bb(:,2)/vox(2))*vox(2);
%bb(:,3) = round(bb(:,3)/vox(3))*vox(3);

M   = prm.VG(1).mat;
vxg = sqrt(sum(M(1:3,1:3).^2));
if det(M(1:3,1:3))<0, vxg(1) = -vxg(1); end;
ogn = M\[0 0 0 1]';
ogn = ogn(1:3)';

% Convert range into range of voxels within template image
x   = (bb(1,1):vox(1):bb(2,1))/vxg(1) + ogn(1);
y   = (bb(1,2):vox(2):bb(2,2))/vxg(2) + ogn(2);
z   = (bb(1,3):vox(3):bb(2,3))/vxg(3) + ogn(3);

og  = -vxg.*ogn;

% Again, chose whether to round to closest voxel.
%of  = -vox.*(round(-bb(1,:)./vox)+1);
of = bb(1,:)-vox;

M1  = [vxg(1) 0 0 og(1) ; 0 vxg(2) 0 og(2) ; 0 0 vxg(3) og(3) ; 0 0 0 1];
M2  = [vox(1) 0 0 of(1) ; 0 vox(2) 0 of(2) ; 0 0 vox(3) of(3) ; 0 0 0 1];
mat = prm.VG(1).mat*inv(M1)*M2;

LEFTHANDED = true;
if (LEFTHANDED && det(mat(1:3,1:3))>0) || (~LEFTHANDED && det(mat(1:3,1:3))<0),
    Flp = [-1 0 0 (length(x)+1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
    mat = mat*Flp;
    x   = flipud(x(:))';
end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [bb,vx] = bbvox_from_V(V)
vx = sqrt(sum(V.mat(1:3,1:3).^2));
if det(V.mat(1:3,1:3))<0, vx(1) = -vx(1); end;

o  = V.mat\[0 0 0 1]';
o  = o(1:3)';
bb = [-vx.*(o-1) ; vx.*(V.dim(1:3)-o)];
return;
%_______________________________________________________________________

%_______________________________________________________________________
function VO = write_dets(P,bb,vox)
if nargin==1,
    job = P;
    P   = job.P;
    bb  = job.bb;
    vox = job.vox;
end;

spm_progress_bar('Init',numel(P),'Writing','volumes completed');

for i=1:numel(V),
    prm = load(deblank(P{i}));
    [x,y,z,mat] = get_xyzmat(prm,bb,vox);
    Tr = prm.Tr;
    BX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1);
    BY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1);
    BZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1);
    DX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1,'diff');
    DY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1,'diff');
    DZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1,'diff');

    [pth,nam,ext,nm] = spm_fileparts(P{i});
    VO = struct('fname',fullfile(pth,['jy_' nam ext nm]),...
        'dim',[numel(x),numel(y),numel(z)],...
        'dt',[spm_type('float32') spm_platform('bigend')],...
        'pinfo',[1 0 0]',...
        'mat',mat,...
        'n',1,...
        'descrip','Jacobian determinants');
    VO     = spm_create_vol(VO);
    detAff = det(prm.VF(1).mat*prm.Affine/prm.VG(1).mat);
    Dat    = single(0);
    Dat(VO.dim(1),VO.dim(2),VO.dim(3)) = 0;

    for j=1:length(z),   % Cycle over planes
        % Nonlinear deformations
        tx = get_2Dtrans(Tr(:,:,:,1),BZ,j);
        ty = get_2Dtrans(Tr(:,:,:,2),BZ,j);
        tz = get_2Dtrans(Tr(:,:,:,3),BZ,j);

        %----------------------------------------------------------------------------
        j11 = DX*tx*BY' + 1; j12 = BX*tx*DY';     j13 = BX*get_2Dtrans(Tr(:,:,:,1),DZ,j)*BY';
        j21 = DX*ty*BY';     j22 = BX*ty*DY' + 1; j23 = BX*get_2Dtrans(Tr(:,:,:,2),DZ,j)*BY';
        j31 = DX*tz*BY';     j32 = BX*tz*DY';     j33 = BX*get_2Dtrans(Tr(:,:,:,3),DZ,j)*BY' + 1;
        % The determinant of the Jacobian reflects relative volume changes.
        %------------------------------------------------------------------
        dat       = (j11.*(j22.*j33-j23.*j32) - j21.*(j12.*j33-j13.*j32) + j31.*(j12.*j23-j13.*j22)) * detAff;
        Dat(:,:,j) = single(dat);
        if numel(P)<5, spm_progress_bar('Set',i-1+j/length(z)); end;
    end;
    VO = spm_write_vol(VO,Dat);
    spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');
return;
%_______________________________________________________________________


function models=train_svms(imnames, feat_dir, region_meta_info, cachedir, run_name, varargin)
%random seed
rng(3);

%input parser to allow additional inputs
p = inputParser;
addOptional(p, 'to_train', [1:20], @isnumeric);
addOptional(p, 'svm_C', 10^-3, @isscalar);
addOptional(p,'bias_mult',10,@isscalar);
addOptional(p, 'pos_loss_weight', 2, @isscalar);
addOptional(p, 'do_latent', 1, @isscalar);
addOptional(p, 'pos_ov_thresh', 0.5, @isscalar);
addOptional(p, 'neg_ov_thresh', 0.2, @isscalar);
addOptional(p, 'hard_neg_thresh', -1.0001, @isscalar);
addOptional(p, 'evict_thresh', -1.2, @isscalar);
addOptional(p, 'retrain_limit', 2000, @isscalar);
parse(p, varargin{:});
opts=p.Results;

to_train=opts.to_train;

%if cachedir does not exist, create it
if(~exist(cachedir, 'file')) mkdir(cachedir); end
cachedir=fullfile(cachedir, run_name);
if(~exist(cachedir, 'file')) mkdir(cachedir); end

%collect feature stats
fprintf('Collecting feature stats...\n');
stats=collect_feature_stats(imnames, feat_dir, cachedir, region_meta_info);

%initialize caches
for i=1:numel(to_train)
    caches(i) = init_cache;
    models(i).categid=to_train(i);
    models(i).w=[];
    models(i).b=[];
    models(i).meannrm=stats.meannrm;
end

%load GT positives
pos_file=fullfile(cachedir, 'pos.mat');
if(exist(pos_file, 'file'))
    tmp=load(pos_file);
    X_pos=tmp.X_pos;
else
    X_pos=get_positives(imnames, feat_dir, to_train, region_meta_info, false, models, opts.pos_ov_thresh);
    save(pos_file, 'X_pos');
end
for i=1:numel(to_train)
    caches(i).X_pos=X_pos{i};
    fprintf('Category %d has %d positives\n', to_train(i), size(X_pos{i},1));
end

%train
first_time=true;
force_update=false;
max_outer_iter=opts.do_latent+1;
for iter=1:max_outer_iter
    if(iter~=max_outer_iter)
        %subsample negatives
        negidx=randperm(numel(imnames));
        negidx=negidx(1:min(1000, numel(negidx)));
    else
        negidx=1:numel(imnames);
    end
    if(iter>1)
        %latent update
        X_pos=get_positives(imnames, feat_dir, to_train, region_meta_info, true, models, opts.pos_ov_thresh);
        for i=1:numel(to_train)
            caches(i).X_pos=X_pos{i};
            fprintf('Latent update: Category %d has %d positives\n', to_train(i), size(X_pos{i},1));

        end
        force_update=true;
    end
    %for every negative
    for i=1:numel(negidx)
        fprintf('Hard negatives iter: %d, image number: %d/%d\n',iter, i, numel(negidx));
        [X_neg, curr_keys]=get_hard_negatives(negidx(i), imnames, feat_dir, to_train,...
                             region_meta_info, models, {caches.keys}, opts.neg_ov_thresh, opts.hard_neg_thresh, first_time);
       very_last_iter=(iter==max_outer_iter) && (i==numel(negidx));
       for j=1:numel(to_train)
            caches(j).X_neg=cat(1, caches(j).X_neg,X_neg{j});
            caches(j).keys=cat(1, caches(j).keys, curr_keys{j});
            caches(j).num_added=caches(j).num_added+size(X_neg{j},1);
            if(caches(j).num_added>opts.retrain_limit || first_time ||very_last_iter || force_update)
                %time to retrain
                fprintf('Retraining category %d with %d positives and %d negatives\n', to_train(j), size(caches(j).X_pos,1), size(caches(j).X_neg,1));
                models(j)=update_model(models(j), caches(j), opts.svm_C, opts.bias_mult, opts.pos_loss_weight);
                [reg, posloss, negloss, neg_scores]=compute_obj_val(models(j),...
                                                   caches(j), opts.svm_C, opts.pos_loss_weight);
                caches(j).pos_loss(end+1)=posloss;
                caches(j).neg_loss(end+1)=negloss;
                caches(j).reg_loss(end+1)=reg;
                caches(j).tot_loss(end+1)=posloss+negloss+reg;
                for t=1:numel(caches(j).tot_loss)
                    fprintf('%d: %f = %f + %f + %f\n',t,caches(j).tot_loss(t),...
                    caches(j).reg_loss(t),caches(j).pos_loss(t), caches(j).neg_loss(t));
                end

                %evict
                keep=neg_scores>opts.evict_thresh;
                caches(j).X_neg=caches(j).X_neg(keep,:);
                caches(j).keys=caches(j).keys(keep,:);
                fprintf('Kept: %d negatives \n', sum(keep));
                caches(j).num_added=0;
            end

        end
        first_time=false;
        force_update=false;
    end
end

%--------------------------------
%   Feature stats
%--------------------------------
function stats = collect_feature_stats(imnames, feat_dir, cachedir, region_meta_info)
stat_file=fullfile(cachedir, 'stats.mat');
if(exist(stat_file, 'file'))
    tmp=load(stat_file);
    stats=tmp.stats; clear tmp;
else

    num_images=100;
    regs_per_img=100;
    ri=randperm(numel(imnames));
    ri=ri(1:num_images);
    totalnrm=0;
    totalcnt=0;
    for i=1:num_images
        fprintf('%d/%d\n', i, num_images);

        %load random image
        tmp=load(fullfile(feat_dir, [imnames{ri(i)} '.mat']));

        %pick random regions
        num=min(regs_per_img, region_meta_info.num_regs(ri(i)));
        ri2=randperm(region_meta_info.num_regs(ri(i)));
        ri2=ri2(1:num);

        %get norm
        nrm=sqrt(sum(tmp.feats(ri2,:).^2,2));
        totalnrm=totalnrm+sum(nrm);
        totalcnt=totalcnt+num;
    end
    stats.meannrm=totalnrm./totalcnt;
    save(stat_file, 'stats');
end
fprintf('\n\nMean norm=%f\n\n',stats.meannrm);

%------------------------------
%   Init cache
%------------------------------
function cache=init_cache
cache.X_pos = single([]);
cache.X_neg = single([]);
cache.keys = [];
cache.num_added = 0;
cache.retrain_limit = 2000;
cache.evict_thresh = -1.2;
cache.hard_thresh = -1.0001;
cache.pos_loss = [];
cache.neg_loss = [];
cache.reg_loss = [];
cache.tot_loss = [];

%------------------------------
%   Get positives from image
%------------------------------

function X_pos=get_positives(imnames, featdir, to_train, region_meta_info, islatent, models, pos_ov_thresh)
for k=1:numel(to_train)
    X_pos{k}=single([]);
end

for i=1:numel(imnames)
    fprintf('Positives: %d/%d\n',i, numel(imnames));
    tmp=load(fullfile(featdir, [imnames{i} '.mat']));
    feats=xform_feat(tmp.feats, models(1).meannrm);

    %for every category
    for j=1:numel(to_train)
        sel=[];
        if(islatent)
            %latent update: consider all features greater than thresh
            %find gt for this category and get corresponding overlaps
            gtidx=find(region_meta_info.gt{i}==to_train(j));
            nongt=find(region_meta_info.gt{i}==0);

            ovs=region_meta_info.overlaps{i}(gtidx,:);
            overlapping=find(max(ovs,[],1)>=pos_ov_thresh);
            ovs=ovs(:,overlapping);

            scores=feats*models(j).w+models(j).b;

            %for every gt, find the highest scoring gt that overlaps by >70
            for k=1:numel(gtidx)
                ovs_this=ovs(k,:);
                idx=find(ovs_this>pos_ov_thresh);
                if(isempty(idx)) continue; end
                [m, argmax]=max(scores(nongt(overlapping(idx))));
                sel=[sel nongt(overlapping(idx(argmax)))];
            end
        else
            sel=find(region_meta_info.gt{i}==to_train(j));
        end
        X_pos{j}=cat(1, X_pos{j},feats(sel,:));
    end
end

%------------------------------
%   Hard negatives
%------------------------------

function [X_neg, curr_keys]=get_hard_negatives(i, imnames, featdir, to_train, region_meta_info, models, keys, neg_ov_thresh, score_thresh, first_time)
tmp=load(fullfile(featdir, [imnames{i} '.mat']));
feats=xform_feat(tmp.feats, models(1).meannrm);

if(~first_time)
    %concatenate all models
    w=cat(2,models(:).w);
    b=cat(2,models(:).b);

    %compute all scores
    scores=bsxfun(@plus,feats*w, b);
else
    scores=(score_thresh+1)*ones(size(feats,1), numel(models));
end

%for every category
for j=1:numel(to_train)
    %get the groundtruths that don't belong to this category
    wronggtidx=find(region_meta_info.gt{i}~=to_train(j) & region_meta_info.gt{i}~=0);
    %get the non-groundtruths that overlap by less than threshold

    nongt=find(region_meta_info.gt{i}==0);
    gtidx=find(region_meta_info.gt{i}==to_train(j));
    if(~isempty(gtidx))
        idx=find(max(region_meta_info.overlaps{i}(gtidx,:),[],1)<neg_ov_thresh);
        nongtidx=nongt(idx);
    else
        nongtidx=nongt;
    end
    allidx=[wronggtidx(:); nongtidx(:)];
    allidx=allidx(scores(allidx,j)>=score_thresh);

    %get the keys for each of these
    tmp_keys=[i*ones(numel(allidx),1) allidx];

    %check that duplicate keys are not added
    keep=~ismember(tmp_keys, keys{j}, 'rows');

    %add these feats
    X_neg{j}=feats(allidx(keep),:);
    curr_keys{j}=tmp_keys(keep,:);
end

%------------------------------
%   Update model
%------------------------------
function model=update_model(model, cache, svm_c, bias_mult, pos_weight)
opts = sprintf('-w1 %.10f -c %.10f -s 3 -B %.10f', ...
                      pos_weight, svm_c, bias_mult);

X=cat(2, cache.X_pos', cache.X_neg');
y=[ones(size(cache.X_pos,1),1); zeros(size(cache.X_neg,1),1)];
fprintf('calling liblinear with:%s\n',opts);
llm = liblinear_train(y, sparse(double(X)), opts, 'col');
model.w = single(llm.w(1:end-1)');
model.b = single(llm.w(end)*bias_mult);

%---------------------------
% Compute objective value
%---------------------------
function [reg, posloss, negloss, neg_scores]=compute_obj_val(model,cache, svm_c, pos_weight)
pos_scores=cache.X_pos*model.w+model.b;
neg_scores=cache.X_neg*model.w+model.b;
posloss=sum(max(1-pos_scores,0))*svm_c*pos_weight;
negloss=sum(max(1+neg_scores,0))*svm_c;
reg=0.5*model.w'*model.w;

%---------------------------
% xform_feat
%---------------------------
function feats=xform_feat(feats, nrm)
feats=feats*20./nrm;


text:
function [hbar, htxt] = barsegstats(segnames,segstats)
%  [hbar htxt] = barsegstats(segnames,segstats)
%
% [segnames segindex segstats] = load_segstats('aseg.stats','bert');
% [hbar htxt] = barsegstats(segnames,segstats(:,3))
%
%

%
% barsegstats.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.4 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

mn = min(segstats(:));
mn = 0;
mx = max(segstats(:));
d = mx-mn;
a = mn+.02*d;

nseg = size(segnames,1);
x = 1:nseg;
hbar = bar(segstats,'r');
for n = 1:nseg
  htxt(n) = text(n+.25,a,segnames(n,:));
  set(htxt(n),'rotation',90);
  set(htxt(n),'fontsize',4)
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [deltaC,deltaS]=gravSolidTideOffset(rMoon,rSun,TT1,TT2)
%%GRAVSOLIDTIDEOFFSET Compute the offsets to add the fully normalized
%                     spherical harmonic gravitational coefficients to
%                     handle the effects of the solid Earth tides. Ocean
%                     tides and other effects need to be added
%                     separately  to get the full offset due to tides.
%                     These offsets are for fully normalized spherical
%                     harmonic tide-free models, not zero-tide models.
%
%INPTUS: rMoon A vector from the center of the Earth in the International
%          Terrestrial Reference System (the WGS-84 coordinate system), an
%          Earth-centered, Earth-fixed coordinate system, to the
%          moon. This can be obtained using the readJPLEphem and GCRS2ITRS
%          functions.
%     rSun A vector from the center of the Earth in the International
%          Terrestrial Reference System (the WGS-84 coordinate system), an
%          Earth-centered, Earth-fixed coordinate system, to the
%          Sun. This can be obtained using the readJPLEphem and GCRS2ITRS
%          functions.
%  TT1,TT2 Two parts of a Julian date given in terrestrial time (TT). The
%          units of the date are days. The full date is the sum of both
%          terms. The date is broken into two parts to provide more bits of
%          precision. It does not matter how the date is split.
%
%OUTPUTS: deltaC An array holding the offsets to the fully normalized
%                coefficient terms that are multiplied by cosines in the
%                spherical harmonic expansion of the gravitational
%                potential. if given to a CountingClusterSet class, then
%                C(n+1,m+1) is the coefficient of degree n and order m. The
%                offsets only go up to degree and order 4.
%         deltaS An array holding the coefficient terms that are multiplied
%                by sines in the spherical harmonic expansion. The format
%                of S is the same as that of C.
%
%The algorithm is the one described in Chapter 5 of [1], except values for
%k_{2m}^{(+)} taken from the 2010 IERS conventions in [2] for an elastic
%Earth model are added, because [1] does not provide the necessary
%coefficients.
%
%The model provided in [1] is not as high order a model as that used in
%Section 6 of [2]. However, the second step of the model in[2] is not easy
%to understand.
%
%To test this function, one can use the same values for the Sun and Moon
%that are given in the IERS's implementation, DEHANTTIDEINEL.F for a
%particular date. Those are
% rSun=[137859926952.015;54228127881.4350;23509422341.6960];
% rMoon=[-179996231.920342;-312468450.131567;-169288918.592160];
% [Jul1,Jul2]=Cal2UTC(2009,4,13,0,0,0);
% [TT1,TT2]=UTC2TT(Jul1,Jul2);
% [deltaC,deltaS]=gravSolidTideOffset(rMoon,rSun,TT1,TT2);
%
%To get the full offset in the gravitational coefficients due to tides,
%the functions gravSolidTideOffset, gravOceanTideOffset and
%gravPoleTideOffset have to be combined.
%
%REFERENCES:
%[1] S. E. Urban and P. K. Seidelmann, Explanatory Supplement to the
%    Astronomical Almanac, 3rd ed., Eds. Mill Valley, CA: University
%    Science Books, 2013.
%[2] G. Petit and B. Luzum, IERS Conventions (2010), International Earth
%    Rotation and Reference Systems Service Std. 36, 2010.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The convention for indexing arrays of constant parameters is that 1 refers
%to a value with respect to lunar parameters and 2 refers to a value with
%respect to solar parameters.

%The equatorial radius of the Earth
Re=Constants.EarthEqRadius;
GMRat(1)=Constants.MoonEarthMassRatio;
GMRat(2)=Constants.SunEarthMassRatio;

%The corrections are only up to the fourth order coefficients.
M=4;
totalNumCoeffs=(M+1)*(M+2)/2;
emptyData=zeros(totalNumCoeffs,1);
deltaC=CountingClusterSet(emptyData);
deltaS=CountingClusterSet(emptyData);

%Obtain the spherical coordinates of the moon.
pointSpherMoon=Cart2Sphere(rMoon);
rj(1)=pointSpherMoon(1);
Lambdaj(1)=pointSpherMoon(2);
Phij(1)=pointSpherMoon(3);

%Obtain the spherical coordinates of the sun.
pointSpherSun=Cart2Sphere(rSun);
rj(2)=pointSpherSun(1);
Lambdaj(2)=pointSpherSun(2);
Phij(2)=pointSpherSun(3);

%%%%These are the Step 1 corrections

%The nomial value from Table 5.5
k2=0.3;

%Implementing Equations 5.74, 5.75, and 5.76
sumVal20=0;
sumVal21=0;
sumVal22=0;
for n=1:2
    P=legendre(2,sin(Phij(n)));
    P20=P(1);
    P21=P(2);
    P22=P(3);

    sumVal20=sumVal20+GMRat(n)*(Re/rj(n))^3*P20;
    sumVal21=sumVal21+GMRat(n)*(Re/rj(n))^3*P21*exp(-1i*Lambdaj(n));
    sumVal22=sumVal22+GMRat(n)*(Re/rj(n))^3*P22*exp(-1i*2*Lambdaj(n));
end
deltaC(2+1,0+1)=(1/sqrt(5))*k2*sumVal20;

deltaC(2+1,1+1)=(1/3)*sqrt(3/5)*k2*real(sumVal21);
deltaS(2+1,1+1)=-(1/3)*sqrt(3/5)*k2*imag(sumVal21);

deltaC(2+1,2+1)=(1/12)*sqrt(12/5)*k2*real(sumVal22);
deltaS(2+1,2+1)=-(1/12)*sqrt(12/5)*k2*imag(sumVal22);

%In the book, no values for k_{2m}^{(+)} were given, so the values here are
%taken from the elastic Earth model in Table 6.3 of the IERS 2010
%conventions.
knmPlus=[-0.00087;-0.00079;-0.00057];
%This loop implements Equation 5.77
for m=0:2

    sumVal=0;
    for n=1:2
        P=legendre(2,sin(Phij(n)));
        P2m=P(m+1);

        sumVal=sumVal+GMRat(n)*(Re/rj(n))^3*P2m*exp(-1i*m*Lambdaj(n));
    end

    deltaC(4+1,m+1)=(1/5)*knmPlus(m+1)*real(sumVal);
    deltaS(4+1,m+1)=-(1/5)*knmPlus(m+1)*imag(sumVal);
end

%%%%Now for the step 2 (frequency-dependent) corrections.
%Convert the terrestrial time to Greenwich mean sidereal time in radians.
GMST=TT2GMST(TT1,TT2);

%These parameters are taken from Table 5.5. The columns are
%Doodson parameters of tau, s, h, p, N', P1 and Am*delta_kx*Hs
m1TideParam=[1, -1,  0,  0,  0,  0,  -16.4e-12;
             1,  1, -2,  0,  0,  0,  -49.6e-12;
             1,  1,  0,  0, -1,  0,   -9.4e-12;
             1,  1,  0,  0,  0,  0,  507.4e-12;
             1,  1,  0,  0,  1,  0,   73.5e-12;
             1,  1,  1,  0,  0, -1,  -15.2e-12];
m2TideParam=[2,  0,  0,  0,  0,  0,   39.5e-12;
             2,  2,  -2, 0,  0,  0,   18.4e-12];

%Turn terrestrial time into barycentric dynamical time TDB.
[TDB1,TDB2]=TT2TDB(TT1,TT2);

%Get the Delaunay variables in radians
vec=DelaunayVar(TDB1,TDB2);
l=vec(1);
lp=vec(2);
F=vec(3);
D=vec(4);
Omega=vec(5);

s=F+Omega;
beta(2)=s;%s: Moon's mean longitude
beta(1)=GMST+pi-s;%tau: time angle in lunar days reckoned from lower transit
beta(3)=s-D;%h: Sun's mean longitude
beta(4)=s-l;%p: longitude of the Moon's mean perigee
beta(5)=-Omega;%N':negative longitude of the Moon's mean node
beta(6)=s-D-lp;%pl: Longitude of the Sun's mean perigee.

%Diurnal tide terms affect the n=2, m=1 coefficients.
n=2;
m=1;
sumVal=0;
for curThetaS=1:size(m1TideParam,1)
    thetaS=sum(beta.*m1TideParam(curThetaS,1:6));

    sumVal=sumVal+m1TideParam(curThetaS,7)*(-1i)*exp(1i*thetaS);
end
deltaC(n+1,m+1)=deltaC(n+1,m+1)+real(sumVal);
deltaS(n+1,m+1)=deltaS(n+1,m+1)-imag(sumVal);

%Semi-durnal tides affects the n=2, m=2 coefficients.
n=2;
m=2;

sumVal=0;
for curThetaS=1:size(m2TideParam,1)
    thetaS=sum(beta.*m2TideParam(curThetaS,1:6));

    sumVal=sumVal+m2TideParam(curThetaS,7)*(-1i)*exp(1i*thetaS);
end
deltaC(n+1,m+1)=deltaC(n+1,m+1)+real(sumVal);
deltaS(n+1,m+1)=deltaS(n+1,m+1)-imag(sumVal);

deltaC=deltaC.clusterEls;
deltaS=deltaS.clusterEls;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function DM = ft_preproc_online_downsample_init(factor)

% function DM = ft_preproc_online_downsample_init(factor)
%
% Initialize an downsampling model with the given factor.

% Copyright (C) 2010, Stefan Klanke
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if factor < 1 or factor~=round(factor)
	error('Argument ''factor'' must be a positive integer number');
end

DM.factor  = factor;
DM.numSkip = 0;


text:
function varargout= bbci_acquire_randomSignals(varargin)
%BBCI_ACQUIRE_RANDOMSIGNALS - Generate random signals
%
%Synopsis:
%  STATE= bbci_acquire_randomSignal('init', <PARAM>)
%  [CNTX, MRKTIME, MRKDESC, STATE]= bbci_acquire_randomSignals(STATE)
%  bbci_acquire_randomSignals('close', STATE)
%
%Arguments:
%  PARAM - Optional arguments:
%    'marker_mode': [CHAR, default 'global'] specifies how markers should be
%             read:
%                '' means no markers,
%                'global' means checking the global variable ACQ_MARKER, and
%                'pyff_udp' means receiving markers by UDP in the pyff format
%    'fs':    [DOUBLE, default 100] sampling rate
%    'clab':  {CELL of CHAR} channel labels,
%             default {'F3','Fz','F4','C3','Cz','C4','P3','Pz','P4'}.
%    'amp': [DOUBLE, default 30] amplitude of the signals (factor for 'randn')
%    'realtime' - [DOUBLE] For value 1, this function only returns the next
%          block of data, if the time corresponding to OPT.blocksize has
%          elaped since the last delivery of data. Prior to that, empty
%          output is returned (just as the 'true' online acquire function
%          would do while no new data is available).
%          For value 0 (default), this function returns one block of data
%          at each call. Values above 1 result in speeded-up
%          realtime. E.g., for 2 the next block of data is returned, if
%          1/2 of the time corresponding to OPT.blocksize has elapsed since
%          the last delivery of data, which amounts to a speed-up factor of 2.
%
%Output:
%  STATE - Structure characterizing the incoming signals; fields:
%     'fs', 'clab', and intern stuff
%  CNTX - 'acquired' signals [Time x Channels]
%  The following variables hold the markers that have been 'acquired' within
%  the current block (if any).
%  MRKTIME - DOUBLE: [1 nMarkers] position [msec] within data block.
%  MRKDESC - DOUBLE: [1 nMarkers] corresponding marker values

% 02-2012 Benjamin Blankertz

% This is just for opt.marker_mode='global':
global ACQ_MARKER

if isequal(varargin{1}, 'init'),
  state= opt_proplistToStruct(varargin{2:end});
  default_clab= {'F3','Fz','F4','C3','Cz','C4','P3','Pz','P4'};
  props= {'fs'            100            '!DOUBLE[1]'
          'clab'          default_clab   'CELL{CHAR}'
          'blocksize'     40             '!DOUBLE[1]'
          'amplitude'     30             '!DOUBLE[1]'
          'marker_mode'   'global'       'CHAR(global pyff_udp)'
          'realtime'      0              '!DOUBLE[1]'
         };
  state= opt_setDefaults(state, props, 1);
  state.nChannels= length(state.clab);
  state.blocksize_sa= ceil(state.blocksize*state.fs/1000);
  state.blocksize= state.blocksize_sa*1000/state.fs;
  state.nsamples= 0;
  if state.realtime==0,
    state.realtime= inf;
  end
  state.start_time= tic;
  %  switch(state.marker_mode),
  %   case 'pyff_udp',
  %    state.socket= open_udp_socket();
  %  end
  output= {state};
elseif isequal(varargin{1}, 'close'),
  if length(varargin)>1,
    state= varargin{2};
    switch(state.marker_mode),
     case 'pyff_udp',
      close_udp_socket(state.socket);
    end
  end
  return
elseif length(varargin)~=1,
  error('Except for INIT/CLOSE case, only one input argument expected');
else
  if isstruct(varargin{1}),
    state= varargin{1};
    time_running= toc(state.start_time);
    if time_running < state.nsamples/state.fs/state.realtime,
      output= {[], [], [], state};
      varargout= output(1:nargout);
      return;
    end
    cntx= state.amplitude*randn(state.blocksize_sa, state.nChannels);
    state.nsamples= state.nsamples + state.blocksize_sa;
    switch(state.marker_mode),
     case '',
      mrkTime= [];
      mrkDesc= [];
     case 'pyff_udp',
      packet= receive_udp(sock);
      if ~isempty(packet),
        % we don't know the marker position within the block -> set randomly
        mrkTime= ceil(state.blocksize_sa*rand)*1000/state.fs;
        mrkDesc= str2int(packet);  % -> check format
      else
        mrkTime= [];
        mrkDesc= [];
      end
     case 'global',
      if ~isempty(ACQ_MARKER),
        mrkTime= ceil(state.blocksize_sa*rand)*1000/state.fs;
        mrkDesc= ACQ_MARKER;
        ACQ_MARKER= [];
      else
        mrkTime= [];
        mrkDesc= [];
      end
     otherwise
      error('unknown marker_mode: %s', state.marker_mode);
    end
    output= {cntx, mrkTime, mrkDesc, state};
  else
    error('unrecognized input argument');
  end
end
varargout= output(1:nargout);

author: bbci


text:
% AR 3/18/2018
function d = dell2ell_pair(trx,fly1,fly2,istry)

nsamples = 20;

%initialize
d = nan(1,trx(fly1).nframes);

%get start and end frames of overlap
t0 = max(trx(fly1).firstframe,trx(fly2).firstframe);
t1 = min(trx(fly1).endframe,trx(fly2).endframe);

% no overlap
if t1 < t0
    return;
end

% ellipse 1

x_mm1 = trx(fly1).x_mm;
y_mm1 = trx(fly1).y_mm;
a_mm1 = trx(fly1).a_mm;
b_mm1 = trx(fly1).b_mm;
theta_mm1 = trx(fly1).theta_mm;

% ellipse 2

x_mm2 = trx(fly2).x_mm;
y_mm2 = trx(fly2).y_mm;
a_mm2 = trx(fly2).a_mm;
b_mm2 = trx(fly2).b_mm;
theta_mm2 = trx(fly2).theta_mm;

off1 = trx(fly1).off;
off2 = trx(fly2).off;

if nargin < 4
    tstry = t0:t1;
else
    % TO DO: doesn't work for two flies?
    tstry = istry(:)' - off1;
end

for t = tstry
    i = t + off1;
    j = t + off2;
    d(i) = ellipse2ellipsedist_hack(x_mm1(i),y_mm1(i), ...
        2*a_mm1(i),2*b_mm1(i),theta_mm1(i), ...
        x_mm2(j),y_mm2(j), ...
        2*a_mm2(j),2*b_mm2(j),theta_mm2(j), ...
        nsamples);
end

author: kristinbranson

repo: JAABA


text:
function [vol, M, mr_parms, volsz] = load_mgh(fname,slices,frames,headeronly)
% [vol, M, mr_parms, volsz] = load_mgh(fname,<slices>,<frames>,<headeronly>)
%
% fname - path of the mgh file
%
% slices - list of one-based slice numbers to load. All
%   slices are loaded if slices is not specified, or
%   if slices is empty, or if slices(1) <= 0.
%
% frames - list of one-based frame numbers to load. All
%   frames are loaded if frames is not specified, or
%   if frames is empty, or if frames(1) <= 0.
%
% M is the 4x4 vox2ras transform such that
% y(i1,i2,i3), xyz1 = M*[i1 i2 i3 1] where the
% indices are 0-based. If the input has multiple frames,
% only the first frame is read.
%
% mr_parms = [tr flipangle te ti fov]
%
% volsz = size(vol). Helpful when using headeronly as vol is [].
%
% See also: save_mgh, vox2ras_0to1
%

%
% load_mgh.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: greve $
%    $Date: 2016/01/19 21:18:27 $
%    $Revision: 1.22 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

vol = [];
M = [];
mr_parms = [];
volsz = [];

if(nargin < 1 | nargin > 4)
  msg = 'USAGE: [vol M] = load_mgh(fname,<slices>,<frames>,<headeronly>)';
  fprintf('%s',msg);
  return;
end

% unzip if it is compressed
if (strcmpi(fname((strlen(fname)-3):strlen(fname)), '.MGZ') | ...
		strcmpi(fname((strlen(fname)-3):strlen(fname)), '.GZ'))
  rand('state', sum(100*clock));
  gzipped =  round(rand(1)*10000000 + ...
		   sum(int16(fname))) + round(cputime);
  %ind = findstr(fname, '.');
  %new_fname = sprintf('/tmp/tmp.load_mgh.%d.mgh', gzipped);
  new_fname = sprintf('%s.load_mgh.%d.mgh', tempname,gzipped);

  if(strcmp(computer,'MAC') || strcmp(computer,'MACI') || ismac)
    unix(sprintf('gunzip -c %s > %s', fname, new_fname)) ;
  else
    unix(sprintf('zcat %s > %s', fname, new_fname)) ;
  end
  fname = new_fname ;
else
  gzipped = -1 ;
end

if(exist('slices')~=1) slices = []; end
if(isempty(slices)) slices = 0; end
if(slices(1) <= 0) slices = 0; end

if(exist('frames')~=1) frames = []; end
if(isempty(frames)) frames = 0; end
if(frames(1) <= 0) frames = 0; end

if(exist('headeronly')~=1) headeronly = 0; end

fid    = fopen(fname, 'rb', 'b') ;
if(fid == -1)
  fprintf('ERROR: could not open %s for reading\n',fname);
  return;
end
v       = fread(fid, 1, 'int') ;
if(isempty(v))
  fprintf('ERROR: problem reading fname\n');
  if(gzipped >=0) unix(sprintf('rm -f %s', fname)); end
end
ndim1   = fread(fid, 1, 'int') ;
ndim2   = fread(fid, 1, 'int') ;
ndim3   = fread(fid, 1, 'int') ;
nframes = fread(fid, 1, 'int') ;
type    = fread(fid, 1, 'int') ;
dof     = fread(fid, 1, 'int') ;

if(slices(1) > 0)
  ind = find(slices > ndim3);
  if(~isempty(ind))
    fprintf('ERROR: load_mgh: some slices exceed nslices\n');
    return;
  end
end

if(frames(1) > 0)
  ind = find(frames > nframes);
  if(~isempty(ind))
    fprintf('ERROR: load_mgh: some frames exceed nframes\n');
    return;
  end
end

UNUSED_SPACE_SIZE= 256;
USED_SPACE_SIZE = (3*4+4*3*4);  % space for ras transform

unused_space_size = UNUSED_SPACE_SIZE-2 ;
ras_good_flag = fread(fid, 1, 'short') ;
if (ras_good_flag)
  delta  = fread(fid, 3, 'float32') ;
  Mdc    = fread(fid, 9, 'float32') ;
  Mdc    = reshape(Mdc,[3 3]);
  Pxyz_c = fread(fid, 3, 'float32') ;

  D = diag(delta);

  Pcrs_c = [ndim1/2 ndim2/2 ndim3/2]'; % Should this be kept?

  Pxyz_0 = Pxyz_c - Mdc*D*Pcrs_c;

  M = [Mdc*D Pxyz_0;  ...
	0 0 0 1];
  ras_xform = [Mdc Pxyz_c; ...
	0 0 0 1];
  unused_space_size = unused_space_size - USED_SPACE_SIZE ;
end

fseek(fid, unused_space_size, 'cof') ;
nv = ndim1 * ndim2 * ndim3 * nframes;
volsz = [ndim1 ndim2 ndim3 nframes];

MRI_UCHAR =  0 ;
MRI_INT =    1 ;
MRI_LONG =   2 ;
MRI_FLOAT =  3 ;
MRI_SHORT =  4 ;
MRI_BITMAP = 5 ;

% Determine number of bytes per voxel
switch type
 case MRI_FLOAT,
  nbytespervox = 4;
 case MRI_UCHAR,
  nbytespervox = 1;
 case MRI_SHORT,
  nbytespervox = 2;
 case MRI_INT,
  nbytespervox = 4;
end

if(headeronly)
  fseek(fid,nv*nbytespervox,'cof');
  if(~feof(fid))
    [mr_parms count] = fread(fid,4,'float32');
    if(count ~= 4)
      fprintf('WARNING: error reading MR params\n');
    end
  end
  fclose(fid);
  if(gzipped >=0)  unix(sprintf('rm -f %s', fname));  end
  return;
end

%------------------ Read in the entire volume ----------------%
if(slices(1) <= 0 & frames(1) <= 0)
  switch type
   case MRI_FLOAT,
    vol = fread(fid, nv, 'float32') ;
   case MRI_UCHAR,
    vol = fread(fid, nv, 'uchar') ;
   case MRI_SHORT,
    vol = fread(fid, nv, 'short') ;
   case MRI_INT,
    vol = fread(fid, nv, 'int') ;
  end

  if(~feof(fid))
    [mr_parms count] = fread(fid,4,'float32');
    if(count ~= 4)
      fprintf('WARNING: error reading MR params\n');
    end
  end
  fclose(fid) ;
  if(gzipped >=0)  unix(sprintf('rm -f %s', fname));  end

  nread = prod(size(vol));
  if(nread ~= nv)
    fprintf('ERROR: tried to read %d, actually read %d\n',nv,nread);
    vol = [];
    return;
  end
  vol = reshape(vol,[ndim1 ndim2 ndim3 nframes]);

  return;
end

%----- only gets here if a subest of slices/frames are to be loaded ---------%

if(frames(1) <= 0) frames = [1:nframes]; end
if(slices(1) <= 0) slices = [1:ndim3]; end

nvslice = ndim1 * ndim2;
nvvol   = ndim1 * ndim2 * ndim3;
filepos0 = ftell(fid);
vol = zeros(ndim1,ndim2,length(slices),length(frames));
nthframe = 1;
for frame = frames

  nthslice = 1;
  for slice = slices
    filepos = ((frame-1)*nvvol + (slice-1)*nvslice)*nbytespervox + filepos0;
    fseek(fid,filepos,'bof');

    switch type
     case MRI_FLOAT,
      [tmpslice nread]  = fread(fid, nvslice, 'float32') ;
     case MRI_UCHAR,
      [tmpslice nread]  = fread(fid, nvslice, 'uchar') ;
     case MRI_SHORT,
      [tmpslice nread]  = fread(fid, nvslice, 'short') ;
     case MRI_INT,
      [tmpslice nread]  = fread(fid, nvslice, 'int') ;
    end

    if(nread ~= nvslice)
      fprintf('ERROR: load_mgh: reading slice %d, frame %d\n',slice,frame);
      fprintf('  tried to read %d, actually read %d\n',nvslice,nread);
      fclose(fid);
      if(gzipped >=0) unix(sprintf('rm -f %s', fname)); end
      return;
    end

    vol(:,:,nthslice,nthframe) = reshape(tmpslice,[ndim1 ndim2]);
    nthslice = nthslice + 1;
  end

  nthframe = nthframe + 1;
end

% seek to just beyond the last slice/frame %
filepos = (nframes*nvvol)*nbytespervox + filepos0;
fseek(fid,filepos,'bof');

if(~feof(fid))
  [mr_parms count] = fread(fid,5,'float32');
  if(count < 4)
    fprintf('WARNING: error reading MR params\n');
  end
end

fclose(fid) ;
if(gzipped >=0) unix(sprintf('rm -f %s', fname)); end

return;


text:
function nlps_master_ex2(alg)
if nargin < 1
    alg = 'DEFAULT';
end
problem = struct( ...
    'f_fcn',    @(x)f2(x), ...
    'gh_fcn',   @(x)gh2(x), ...
    'hess_fcn', @(x, lam, cost_mult)hess2(x, lam, cost_mult), ...
    'x0',       [1; 1; 0], ...
    'opt',      struct('verbose', 2, 'alg', alg) ...
);
[x, f, exitflag, output, lambda] = nlps_master(problem);
fprintf('\nf = %g   exitflag = %d\n', f, exitflag);
fprintf('\nx = \n');
fprintf('   %g\n', x);
fprintf('\nlambda.ineqnonlin =\n');
fprintf('   %g\n', lambda.ineqnonlin);

function [f, df, d2f] = f2(x)
f = -x(1)*x(2) - x(2)*x(3);
if nargout > 1          %% gradient is required
    df = -[x(2); x(1)+x(3); x(2)];
    if nargout > 2      %% Hessian is required
        d2f = -[0 1 0; 1 0 1; 0 1 0];   %% actually not used since
    end                                 %% 'hess_fcn' is provided
end

function [h, g, dh, dg] = gh2(x)
h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
dh = 2 * [x(1) x(1); -x(2) x(2); x(3) x(3)];
g = []; dg = [];

function Lxx = hess2(x, lam, cost_mult)
if nargin < 3, cost_mult = 1; end   %% allows to be used with 'fmincon'
mu = lam.ineqnonlin;
Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + ...
        [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower