text:
function demo_LocNet_object_detection_pipeline(varargin)
% demo of the object detection pipeline using as a localization model
% either a LocNet model or a CNN-based bounding box regression model
%
% The input arguments are given in the form of Name,Value pair arguments
% and are:
% 'gpu_id': scalar value with gpu id (one-based index) that will be used for
% running the demo. If a non positive value is given then the CPU will be
% used instead. Default value: 1
% 'rec_dir_name': the directory name of the recognition model that will be
% used for running the demo. The possible options (if you download and
% properly place the pre-trained models) are 'VGG16_Reduced_MRCNN' for the
% Reduced MR-CNN model and 'VGG16_FastRCNN' for the (re-implemented by us)
% Fast RCNN model. Default value: 'VGG16_Reduced_MRCNN'
% 'loc_dir_name': the directory name of the localization model that will be
% used for running the demo. The possible options (if you download and
% properly place the pre-trained models) are:
% 1) 'VGG16_LocNet_InOut':    LocNet-InOut localization model
% 2) 'VGG16_LocNet_Borders':  LocNet-Borders localization model
% 3) 'VGG16_LocNet_Combined': LocNet-Combined localization model
% 4) 'VGG16_BBoxReg':         CNN-based bounding box regression (baseline)
% Default value: 'VGG16_LocNet_InOut'
% 'box_method': string with the box proposals algorithm that will be used
% in order to generate the set of candidate boxes that will be given as
% input to the object detection pipeline. Currently the possible choises
% are: 'edge_boxes', 'selective_search', or 'dense_boxes_10k'. Default
% value: 'edge_boxes'.
% 'img_path': path to the image that will be used for the demo. Note that
% if you use other than the default image for running the demo you will
% need to first install the box proposal algorithms (read the README.md
% file).
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

ip = inputParser;
ip.addParamValue('gpu_id',      1,                                  @isscalar);
ip.addParamValue('box_method',  'edge_boxes',                       @ischar);
ip.addParamValue('img_path',    './code/examples/images/000084.jpg',@ischar);
ip.addParamValue('rec_dir_name','VGG16_Reduced_MRCNN',              @ischar);
ip.addParamValue('loc_dir_name','VGG16_LocNet_InOut',               @ischar);
ip.parse(varargin{:});
opts = ip.Results;

%******************************* OPTIONS **********************************
% The GPU that will be used for running the demo. gpu_id is a one-based
% index; if a non positive value is given then the CPU will be used instead.
gpu_id = opts.gpu_id;

% Set the directory name of bounding box recognition model; this could be
% the 'VGG16_Reduced_MRCNN' for the Reduced MR-CNN model and
% 'VGG16_FastRCNN' for the (re-implemented by us) Fast RCNN model.
model_rec_dir_name = opts.rec_dir_name;

% Set the directory name of the LocNet localization model or the CNN-based
% bounding box regression model); the possible options (after downloading
% the models) could be:
% 1) 'VGG16_LocNet_InOut':    LocNet-InOut localization model
% 2) 'VGG16_LocNet_Borders':  LocNet-Borders localization model
% 3) 'VGG16_LocNet_Combined': LocNet-Combined localization model
% 4) 'VGG16_BBoxReg':         CNN-based bounding box regression (baseline)
model_loc_dir_name = opts.loc_dir_name;

img_path   = opts.img_path; % path to the image

box_method = opts.box_method; % string with the box proposals algorithm that
% will be used in order to generate the set of candidate boxes that will be
% given as input to the object detection pipeline. Currently the possible
% choises are: 'edge_boxes', 'selective_search', or 'dense_boxes_10k'

% number of iterations that the iterative localization scheme is performed.
num_iterations = 4;

% maximum number of regions/bboxes that can be given in one go in the
% caffe networks such that everything can fit in the GPU memory;
% For a GPU with 6 Gbytes memory:
model_obj_rec_max_rois_num_in_gpu = 200;
model_obj_loc_max_rois_num_in_gpu = 200;
%**************************************************************************

%***************************** LOAD MODELS ********************************
caffe_set_device( gpu_id );
caffe.reset_all();
% Normaly, the loading of the models needs to performed only once before
% processing a bunch of images
fprintf('Loading detection models... '); th = tic;

% set the path of the bounding box recognition model for object detection
full_model_rec_dir  = fullfile(pwd, 'models-exps', model_rec_dir_name); % full path to the model's directory
use_detection_svms  = true;
model_rec_mat_name  = 'detection_model_svm.mat'; % model's matlab filename
full_model_rec_path = fullfile(full_model_rec_dir, model_rec_mat_name); % full path to the model's matlab file
assert(exist(full_model_rec_dir,'dir')>0);
assert(exist(full_model_rec_path,'file')>0);

% Load the recognition moddel
ld = load(full_model_rec_path, 'model');
model_obj_rec = ld.model;
clear ld;
model_obj_rec = load_object_recognition_model_on_caffe(...
    model_obj_rec, use_detection_svms, full_model_rec_dir);

% Set the path of the LocNet localization model that will be used
full_model_loc_dir  = fullfile(pwd, 'models-exps', model_loc_dir_name);
model_loc_mat_name  = 'localization_model.mat'; % model's matlab filename
full_model_loc_path = fullfile(full_model_loc_dir, model_loc_mat_name);  % full path to the model's matlab file
assert(exist(full_model_loc_dir,'dir')>0);
assert(exist(full_model_loc_path,'file')>0);

% Load the LocNet model
ld = load(full_model_loc_path, 'model');
model_obj_loc = ld.model;
clear ld;
model_obj_loc = load_bbox_loc_model_on_caffe(model_obj_loc, full_model_loc_dir);
fprintf(' %.3f sec\n', toc(th));

% maximum number of regions/bboxes that can be given in one go in the
% caffe networks such that everything can fit in the GPU memory
model_obj_rec.max_rois_num_in_gpu = model_obj_rec_max_rois_num_in_gpu;
model_obj_loc.max_rois_num_in_gpu = model_obj_loc_max_rois_num_in_gpu;
%**************************************************************************

%************************** CONFIGURATION *********************************
category_names = model_obj_rec.classes; % a C x 1 cell array with the name
% of the categories that the detection system looks for. C is the number of
% categories.
num_categories = length(category_names);

conf = struct;
conf.thresh = -3 * ones(num_categories,1); % It contains the
% threshold per category that will be used for removing scored boxes with
% low confidence prior to applying the non-max-suppression step.
conf.nms_iou_thrs = 0.3; % IoU threshold for the non-max-suppression step
conf.box_method = box_method; % string with the box proposals algorithm that
% will be used in order to generate the set of candidate boxes that will be
% given as input to the object detection pipeline. Currently the possible
% choises are: 'edge_boxes', 'selective_search', or 'dense_boxes_10k'
conf.num_iterations = num_iterations; % number of iterations that the iterative localization scheme is performed.
conf.thresh_init = -2.5 * ones(num_categories,1); % the threshold per category
% that will be used in order to prune the candidate boxes with low
% confidence only at the first iteration of the object detection pipeline
% (in order to remove the computation cost of subsequent steps). Normally
% this threshold is chosen such as the average number of candidate boxes
% per image and per category to be around 18.
conf.nms_iou_thrs_init = 0.95;
if strcmp(conf.box_method,'dense_boxes_10k'), conf.nms_iou_thrs_init = 0.85; end
% nms_iou_thrs_init: the IoU threshold that will be used during the
% non-max-suppression step that is applied at the first iteration of the
% object detection pipeline in order to remove near duplicate box proposals
%**************************************************************************

%********************* RUN THE DETECTION PIPELINE *************************
img = imread(img_path); % load image
% extract the box proposals that will be given as input to the object
% detection pipeline
fprintf('Extracting bounding box proposals... '); th = tic;
[a,b] = fileparts(img_path);
proposals_files = fullfile(a, [b, '_',box_method, '_box_proposals.mat']);
try
    ld = load(proposals_files);
    bbox_proposals = ld.bbox_proposals;
catch
    bbox_proposals = extract_object_proposals( img, conf );
%     save(proposals_files,'bbox_proposals');
end
fprintf(' %.3f sec\n', toc(th));

% run the LocNet object detection pipeline in the image
fprintf('LocNet object detection pipeline... '); th = tic;
bbox_detections = LocNet_object_detection( img, model_obj_rec, ...
    model_obj_loc, bbox_proposals, conf);
fprintf(' %.3f sec\n', toc(th));

caffe.reset_all(); % free the memory occupied by the caffe models
%**************************************************************************

%************************ VISUALIZE DETECTIONS ****************************
% visualize the bounding box detections.
score_thresh = 0 * ones(num_categories, 1); % the minimum score threshold per
% category for keeping or discarding a detection. For the purposes of this
% demo we set the thresholds to 0 value (which is not however, the
% optimal value).
display_bbox_detections( img, bbox_detections, score_thresh, category_names );
%**************************************************************************
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
% this simple exampe shows the general principles of geodesic toolbox
% Danil Kirsanov, 09/2007

global geodesic_library;
geodesic_library = 'geodesic_debug';      %"release" is faster and "debug" does additional checks
rand('state', 0);                         %comment this statement if you want to produce random mesh every time

N = 300;                                  %number of points in a mesh
[vertices,faces] = create_hedgehog_mesh(N, 0.1);   %create "noisy sphere" mesh; "vertices" contains 3D vertex coordinates; "faces" contains vertex id's for every triangle
%[vertices,faces] = create_flat_triangular_mesh(0.1, 0); N = length(vertices);  %rectangular mesh for sanity check

mesh = geodesic_new_mesh(vertices,faces);         %initilize new mesh
algorithm = geodesic_new_algorithm(mesh, 'exact');      %initialize new geodesic algorithm

vertex_id = 1;                             %create a single source at vertex #1
source_points = {geodesic_create_surface_point('vertex',vertex_id,vertices(vertex_id,:))};

geodesic_propagate(algorithm, source_points);   %propagation stage of the algorithm (the most time-consuming)

vertex_id = N;                              %create a single destination at vertex #N
destination = geodesic_create_surface_point('vertex',vertex_id,vertices(vertex_id,:));
path = geodesic_trace_back(algorithm, destination);     %find a shortest path from source to destination

distances = zeros(N,1);              %find distances to all vertices of the mesh (actual pathes are not computed)

[source_id, distances] = geodesic_distance_and_source(algorithm);     %find distances to all vertices of the mesh; in this example we have a single source, so source_id is always equal to 1

geodesic_delete;                            %delete all meshes and algorithms

%-----------------plotting------------------------
hold off;
colormap('default');
trisurf(faces,vertices(:,1),vertices(:,2),vertices(:,3),distances, 'FaceColor', 'interp', 'EdgeColor', 'k');       %plot the mesh
daspect([1 1 1]);

hold on;
plot3(source_points{1}.x, source_points{1}.y, source_points{1}.z, 'or', 'MarkerSize',3);    %plot sources

plot3(destination.x, destination.y, destination.z, 'ok', 'MarkerSize',3);       %plot destination
[x,y,z] = extract_coordinates_from_path(path);                                  %prepare path data for plotting
h = plot3(x*1.001,y*1.001,z*1.001,'k-','LineWidth',2);    %plot path
legend(h,'geodesic curve');

author: Sable

repo: mcbench-benchmarks


text:
function a = tanh(a)
%TANH         Gradient hyperbolic tangent tanh(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    iaccelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  ax = 1 ./ sqr(cosh(full(a.x(:))));
  a.x = tanh(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function mae = CalMAE(smap, gtImg)
% Code Author: Wangjiang Zhu
% Email: wangjiang88119@gmail.com
% Date: 3/24/2014
if size(smap, 1) ~= size(gtImg, 1) || size(smap, 2) ~= size(gtImg, 2)
    error('Saliency map and gt Image have different sizes!\n');
end

if ~islogical(gtImg)
    gtImg = gtImg(:,:,1) > 128;
end

smap = im2double(smap(:,:,1));
fgPixels = smap(gtImg);
fgErrSum = length(fgPixels) - sum(fgPixels);
bgErrSum = sum(smap(~gtImg));
mae = (fgErrSum + bgErrSum) / numel(gtImg);

repo: Saliency-Evaluation-Toolbox


text:
 function ob = Jop(varargin)
%function ob = Jop(varargin)
%J = Jop({'J', 4, 'L', 4, 'kspace', kspace, 'fov', fov, ...
%     'kn.ktype', 'kaiser', 'kn.kb_alf', kb_alf, 'kn.kb_m', kb_m});
%
% Creates fatrix J for use in operations of the form y = J * x.  Designed
% for use in iterative algorithms for estimating density compensation
% weights for conjugate phase MRI reconstruction.
%
% J is an MxM symmetric banded matrix, with entries corresponding to a
% gridding kernel, C, convolved with itself, evaluated at values
% corresponding to the difference between kspace sample locations, i.e.
%       J(i,l) = C(*)C(v_i - v_l)
% where C is the gridding kernel (usually a Kaiser-Bessel), (*) represents
% the convolution operator, and kspace sample locations are {v_m},
% m = 1,..., M.
%
% in:
%   kspace          kspace sample locations
%   fov             Field of View
%
%  optional inputs:
%   J               neighborhood for interpolation kernel
%   L               samples per integer
%   kn.*            various information to define kernel
%                   e.g. kn.ktype = 'kaiser', kn.kb_alf = 10, kn.kb_m = 0
% out:
%   J   [nd nd]     Fatrix operator
%
% K. Khalsa, Mar. 2006

% default assignments:

arg.J = 4;
arg.L = 4;
arg.kspace = [];
arg.fov = [];
arg.del = [];

arg.kn.ktype = 'kaiser';
arg.kn.kernel = {};     %*** find out what to do with this, see Gn.arg.st
arg.kn.kb_alf = [];     %
arg.kn.kb_m = [];       %

% pair input arguments
if iscell(varargin)
    varargs = varargin{1};
    arg = vararg_pair(arg, varargs);
else
    printf('input needs to be in a cell array.');
    help Jop;
    return;
end

if isempty(arg.kspace)
    error('kspace sample locations are a required input argument');
else
    arg.M = size(arg.kspace, 1);
    arg.dim = [arg.M arg.M];
end

arg.del = 1 ./ (arg.fov * arg.L);

kappa = linspace(-arg.J/2, arg.J/2, arg.J * arg.L + 1)';
%kappa = [-arg.J/2 : 1/arg.L : arg.J/2]';   %equivalently

if (arg.kn.ktype == 'kaiser')

    if xor(isempty(arg.kn.kb_alf), isempty(arg.kn.kb_m))
        printf('For Kaiser-Bessel kernel, both m and alpha are required, or neither');
        return;
    elseif (isempty(arg.kn.kb_alf) && isempty(arg.kn.kb_m))
        % use parameters from Jackson's 1991 gridding paper
        arg.kn.kb_m = 0;
        ww = [1.5 2.0 2.5 3.0 3.5 4.0];
        jack_alf = [6.6875 9.1375 11.5250 13.9086 16.2734 18.5547];
        arg.kn.kb_alf = interp1(ww, jack_alf, arg.J, 'linear', 'extrap');

    end
    % keyboard
    arg.kn.C = kaiser_bessel(kappa, arg.J, arg.kn.kb_alf, ...
        arg.kn.kb_m);
    if size(arg.kspace, 2) == 2
        arg.kn.C = arg.kn.C * arg.kn.C';    % 2D kernel
    end
    arg.kn.c0 = kaiser_bessel_ft(0, arg.J, arg.kn.kb_alf, arg.kn.kb_m);

elseif (arg.kn.ktype == 'function_handle')
    % figure out what to do here
    % see Gn.arg.st.kernel = {[1x1] function_handle}... ??
else
    error('kn.ktype must be kaiser or function_handle');
end

% Build Fatrix object
ob = Fatrix(arg.dim, arg, 'forw', @J_forw, 'back', @J_forw, ...
    'caller', mfilename);
% transpose multiplication = same as fwd multiplication b/c we
% stipulate that all kernels must be real and symmetric

%--------------------------
% multiplication
function y = J_forw(arg, x)

if (size(x,1) ~= arg.dim(2))
    error('dimension mismatch in matrix vector multiplication');
end

% 1D case
if size(arg.kspace, 2) == 1
    k0 = abs(min(arg.kspace / arg.del) - arg.J);
    k = round(arg.kspace / arg.del + k0);
    N = max(k) + arg.J;
    xtmp = full(sum(sparse(1:arg.M, k, x, arg.M, N), 1))';

    CC = conv(arg.kn.C, arg.kn.C);
    Jxtmp = conv(xtmp, CC);
    Jxtmp = arg.del * Jxtmp;
    Jxtmp = Jxtmp / max(CC);    % adjust scaling?
    y = Jxtmp(k + arg.J * arg.L);

    % 2D case
elseif size(arg.kspace, 2) == 2

    k01 = abs(min(arg.kspace(:,1) / arg.del(1)) - arg.J);
    k02 = abs(min(arg.kspace(:,2) / arg.del(2)) - arg.J);
    k0 = [k01, k02];

    k1 = round(arg.kspace(:,1) / arg.del(1) + k0(1));
    k2 = round(arg.kspace(:,2) / arg.del(2) + k0(2));
    N1 = max(k1); N2 = max(k2);

    xtmp = full(sparse(k1, k2, x, N1, N2));

    CC = conv2(arg.kn.C, arg.kn.C);
    Jxtmp = conv2(xtmp, CC, 'same');
    ind = sub2ind(size(Jxtmp), k1,k2);

    y = Jxtmp(ind);
    %  y = y / max(CC(:));  % adjust scaling?
else
    error('only 1D and 2D currently supported')
end
y = y / ((arg.kn.c0)^2); % adjust scaling??
% keyboard

% %plot to see if it's working right 2D
%     figure(11), clf
%     subplot(211)
%     plot3(arg.kspace(:,1), arg.kspace(:,2), x, 'x');
%     title('ks locs before regridding'), xlabel('kx'), ylabel('ky')
%     subplot(212)
%     plot3(k1, k2, x, 'x'), title('ks locs after regridding'), xlabel('kapx')
%
%
%   figure(12), clf
%   subplot(211), imagesc(xtmp), axis square, colorbar, title('xtmp')
%   subplot(212), imagesc(Jx), axis square, colorbar, title('Jx')
% keyboard

%%plot to see if it's working right 1D
%     figure(11), clf
%     subplot(211)
%     stem(arg.kspace, x), title('wi before regridding'), xlabel('kx')
%     subplot(212)
%     stem(1:N, xtmp), title('wi after regridding'), xlabel('kap')

%   figure(12), clf
%     subplot(211), stem(1:length(Jxkap), Jxkap), title('Jxkap')
%     subplot(212), stem(arg.kspace, y), title('y = J * x'), xlabel('kx')
% keyboard

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function  test_softmax_classifier()

    clc;
    clear;
    close all;

    %% Set algorithms
    if 0
        algorithms = gd_solver_list('ALL');
    else
        algorithms = gd_solver_list('LS');
    end

    %% prepare dataset
    if 1
        n_per_class = 100;    % # of samples
        d = 3;      % # of dimensions
        l = 5;      % # of classes
        std = 0.15; % standard deviation

        data = multiclass_data_generator(n_per_class, d, l, std);
        n = length(data.y_train);
        d = d + 1; % adding '1' row for intersect

        % train data
        x_train = [data.x_train; ones(1,n)];
        % transform class label into label logical matrix
        y_train = zeros(l,n);
        for j=1:n
            y_train(data.y_train(j),j) = 1;
        end

        % test data
        x_test = [data.x_test; ones(1,n)];
        % transform class label into label logical matrix
        y_test = zeros(l,n);
        for j=1:n
            y_test(data.y_test(j),j) = 1;
        end

        lambda = 0.0001;
        w_opt = zeros(d*l,1);

    else
        % load real-world data
        data = importdata('../data/mnist/6000_data_0.001.mat');
        x_train = data.x_trn;
        y_train = data.y_trn;
        x_test = data.x_tst;
        y_test= data.y_tst;
        d = size(x_train,1);
        n = length(y_train);
        lambda = data.lambda;

        w_opt = data.w_opt;
        l = data.L;
    end

    % set plot_flag
    if d > 4
        plot_flag = false;  % too high dimension
    else
        plot_flag = true;
    end

    %% define problem definitions
    problem = softmax_regression(x_train, y_train, x_test, y_test, l, lambda);

    % initialize
    w_init = randn(d*l,1);

    w_list = cell(length(algorithms),1);
    info_list = cell(length(algorithms),1);

    %% calculate solution
    if norm(w_opt)
    else
        % calculate solution
        w_opt = problem.calc_solution(1000);
    end
    f_opt = problem.cost(w_opt);
    fprintf('f_opt: %.24e\n', f_opt);

    %% perform algorithms
    for alg_idx=1:length(algorithms)
        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_gnorm = 1e-10;
        options.max_iter = 100;
        options.verbose = true;
        options.f_opt = f_opt;
        options.store_w = true;

        switch algorithms{alg_idx}
            case {'SD-STD'}

                options.step_alg = 'fix';
                options.step_init = 1;
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-SCALE-EXACT'}

                options.sub_mode = 'SCALING';
                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'Newton-STD'}

                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-DAMP'}

                options.sub_mode = 'DAMPED';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-CHOLESKY'}

                options.sub_mode = 'CHOLESKY';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'CG-PRELIM'}

                options.sub_mode = 'PRELIM';
                options.step_alg = 'exact';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-BKT'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'backtracking';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-EXACT'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'exact';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-PRECON-EXACT'}

                options.sub_mode = 'PRECON';
                % diagonal scaling
                options.M = diag(diag(A));
                options.step_alg = 'exact';
                options.beta_alg = 'PR';

                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'NCG-BTK'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'backtracking';
                options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = ncg(problem, options);

            case {'NCG-WOLFE'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'strong_wolfe';
                options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = ncg(problem, options);

            case {'BFGS-H-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-H-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-B-BKT'}

                options.step_alg = 'backtracking';
                options.update_mode = 'B';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-B-EXACT'}

                options.step_alg = 'exact';
                options.update_mode = 'B';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'DAMPED-BFGS-BKT'}

                options.step_alg = 'backtracking';
                options.update_mode = 'Damping';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'DAMPED-BFGS-EXACT'}

                options.step_alg = 'exact';
                options.update_mode = 'Damping';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'L-BFGS-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'BB'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = bb(problem, options);

            case {'SGD'}

                options.batch_size = 1;
                options.step = 0.1 * options.batch_size;
                %options.step_alg = 'decay';
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = sgd(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    fprintf('\n\n');

    %% plot all
    close all;

    % display iter vs cost/gnorm
    display_graph('iter','cost', algorithms, w_list, info_list);
    display_graph('iter','gnorm', algorithms, w_list, info_list);

    % draw convergence sequence
    w_history = cell(1);
    cost_history = cell(1);
    for alg_idx=1:length(algorithms)
        w_history{alg_idx} = info_list{alg_idx}.w;
        cost_history{alg_idx} = info_list{alg_idx}.cost;
    end
    draw_convergence_sequence(problem, w_opt, algorithms, w_history, cost_history);

    % display classification results
    y_pred_list = cell(length(algorithms),1);
    accuracy_list = cell(length(algorithms),1);
    for alg_idx=1:length(algorithms)
        % predict class
        y_pred_list{alg_idx} = problem.prediction(w_list{alg_idx});
        % calculate accuracy
        accuracy_list{alg_idx} = problem.accuracy(y_pred_list{alg_idx});
        fprintf('Classificaiton accuracy: %s: %.4f\n', algorithms{alg_idx}, accuracy_list{alg_idx});
    end

    % convert logial matrix to class label vector
    [~, y_train] = max(y_train, [], 1);
    [~, y_test] = max(y_test, [], 1);
    if plot_flag
        display_classification_result(problem, algorithms, w_list, y_pred_list, accuracy_list, x_train, y_train, x_test, y_test);
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function [D]=logic2levelset(varargin)

% function [D]=logic2levelset(logicInside,voxelSize,voxelSizeResample)
% ------------------------------------------------------------------------
% This function converts the logic image logicInside to a level set image.
% The level set image intensities define the distance of the voxels to the
% boundaries of the logic image.
%
%
% Change log:
%
% ------------------------------------------------------------------------
%%

switch nargin
    case 1
        logicInside=varargin{1};
        voxelSize=[];
        voxelSizeResample=[];
    case 2
        logicInside=varargin{1};
        voxelSize=varargin{2};
        voxelSizeResample=[];
    case 3
        logicInside=varargin{1};
        voxelSize=varargin{2};
        voxelSizeResample=varargin{3};
end

if isempty(voxelSize)
    voxelSize=ones(1,3);
end

if isempty(voxelSizeResample)
    voxelSizeResample=mean(voxelSize);
end

if numel(voxelSizeResample)==1
    voxelSizeResample=voxelSizeResample*ones(1,3);
end

%%

siz=size(logicInside);
%Resample input image isotropically (isotropic voxels)
if ~max(voxelSize-voxelSizeResample)< max(eps(voxelSize))
    resampleOn=1;
    [logicInside]=imageResample(logicInside,voxelSize,voxelSizeResample);
    logicInside=logicInside>0; %Forces logic and fixes potential NaN's
    scaleFactor=voxelSizeResample(1);
else
    resampleOn=0;
    scaleFactor=voxelSize(1);
end
logicInside=logicInside>0; %Force to be a logic (in case resampling altered it)

%Remove interior from logic
try
    logicOn = bwmorph3(logicInside,'remove'); %New in R2018a
catch
    logicOn=logicRemoveInterior(logicInside); %GIBBON alterative
end

%Do distance transform on isotropic image
D = double(bwdist(logicOn,'euclidean')); %Compute distance
D(logicInside)=-D(logicInside); %Negate distance inside
D=D.*scaleFactor; %Scale by voxel size

if resampleOn==1
    sizNew=size(D); %Size of distance data in resampled state
    [J,I,K]=meshgrid(1:1:siz(2),1:1:siz(1),1:1:siz(3)); %Image coordinate grid of original image
    [X,Y,Z]=im2cart(I,J,K,voxelSize); %Spatial coordinate of original image grid
    [I,J,K]=cart2im(X,Y,Z,voxelSizeResample); %Image coordinates of original grid in resampled image
    I=round(I); J=round(J); K=round(K); %Rounding
    IND=reshape(sub2indn(sizNew,[I(:) J(:) K(:)],1),size(I)); %Convert to linear indices
    D=D(IND); %Override D to be data sampled at original image points
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [retval] = freesurfer_fread3(fid)

% freesurfer_fread3 - read a 3 byte integer out of a file
%
% [retval] = freesurfer_fread3(fid)
%
% see also freesurfer_write3, freesurfer_read_surf, freesurfer_write_surf
%

b1 = fread(fid, 1, 'uchar') ;
b2 = fread(fid, 1, 'uchar') ;
b3 = fread(fid, 1, 'uchar') ;
retval = bitshift(b1, 16) + bitshift(b2,8) + b3 ;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function xycoords = latlon2xy(latlon,centerpoint)
    % convert lat,lon list into x,y given a lat-lon centerpoint
    %
    %This function will take coordinates given in latitude and longitude
    %and convert them to and x,y coordinate system with x pointing east, y
    %pointing north.

    %The zero point will be set at the point given in centerpoint, where
    %centerpoint is entered as [lat lon].

    lat = latlon(:,1);
    lon = latlon(:,2);

    %Defining the y and x coordinates of each lat and lon

    y = deg2km(lat - centerpoint(1));

    x = deg2km(lon - centerpoint(2)).*cosd(lat);

    %Outputting a vector with column1 = x coords, column2 = y coords

    xycoords = [x(:) y(:)];
end

author: CelsoReyes


text:
function [ mask ] = sicd_mask( sicd_meta, dim1_range, dim2_range )
%SICD_MASK Create a mask of valid data from SICD metadata
%
% SICD allows a ValidData construct that defines which pixels are valid
% values with a list of vertices of an enclosing polygon.  This functions
% converts that structure into a mask over any portion of the SICD image
% space.
%
% INPUTS:
%   sicd_meta  - required: SICD metadata structure.
%   dim1_range - optional: 1x2 array [first element, last element].  If
%        producing a mask for a subimage, this is the range of image area
%        in the first dimension (what is called "columns" in SICD-- but
%        not necessarily stored in that orientation in the MATLAB reader).
%        Default is entire image range. Can also use an empty array ([]) to
%        specify entire range.
%   dim2_range - optional: 1x2 array [first element, last element].  If
%        producing a mask for a subimage, this is the range of image area
%        in the second (range, called "rows" in SICD) dimension.  Default
%        is entire image range. Can also use an empty array ([]) to specify
%        entire range.
%
% OUTPUTS:
%   mask - Logical array over the area requested
%
% Author: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Default input parameters
if ~exist('row_range','var') || isempty(dim2_range)
    dim2_range = [1 sicd_meta.ImageData.FullImage.NumRows];
end
if ~exist('col_range','var') || isempty(dim1_range)
    dim1_range = [1 sicd_meta.ImageData.FullImage.NumCols];
end

if isfield(sicd_meta.ImageData,'ValidData') && ...
        isfield(sicd_meta.ImageData.ValidData,'Vertex')
    % Poly2mask has some peculiarities with including vertices (which are
    % document in the MATLAB help), so its possible that not all verticies
    % will be included in the resulting mask, but the boundary of the mask
    % will never be more than one pixel off of where it should be.
    mask = poly2mask(...
        [sicd_meta.ImageData.ValidData.Vertex.Row] + 1 ... % ValidData vertices are zero-based
           - dim2_range(1) + 1, ... % Relative to area of interest
        [sicd_meta.ImageData.ValidData.Vertex.Col] + 1 ... % ValidData vertices are zero-based
           - dim1_range(1) + 1, ... % Relative to area of interest
        diff(dim1_range)+1, diff(dim2_range)+1);
else % If valid data is not defined, assume all data is valid.
    mask = true(diff(dim1_range)+1, diff(dim2_range)+1);
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function k = simXsimKernDiagCompute(simKern1, simKern2, t)

% SIMXSIMKERNDIAGCOMPUTE Diagonal of a cross kernel between two SIM kernels.
% FORMAT
% DESC computes diagonal of cross kernel terms between two SIM kernels for
% the multiple output kernel.
% ARG simKern1 : the kernel structure associated with the first SIM
% kernel.
% ARG simKern2 : the kernel structure associated with the second SIM
% kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN k : block of values from the diagonal kernel matrix.
% RETURN sK : normalised diagonal(i.e. unscaled version of K not multiplied
% by sqrt(simKern1.variance), sqrt(simKern2.variance) and
% sqrt(2/simKern.inverseWidth) in the case of the un-normalised kernel).
%
% SEEALSO : multiKernParamInit, multiKernCompute, simKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if size(t, 2) > 1
  error('Input can only have one column');
end

if simKern1.inverseWidth ~= simKern2.inverseWidth
  error('Kernels cannot be cross combined if they have different inverse widths.')
end

if ~isfield(simKern1, 'isNormalised')
    isSim1Normalised = false;
else
    isSim1Normalised = simKern1.isNormalised;
end

if ~isfield(simKern2, 'isNormalised')
    isSim2Normalised = false;
else
    isSim2Normalised = simKern2.isNormalised;
end

% The factor of 2 on the top is because all our derivations are in terms
% of erfs which are defined in terms of exp(-x^2) rather than exp(-0.5x^2).
sigma = sqrt(2/simKern1.inverseWidth);

if ~simKern1.isStationary
    h1 = simXsimComputeDiagH(t, simKern1.decay, simKern2.decay, ...
                     simKern1.delay, simKern2.delay, sigma);
    h2 = simXsimComputeDiagH(t, simKern2.decay, simKern1.decay, ...
                     simKern2.delay, simKern1.delay, sigma);
else
    h1 = simXsimComputeDiagHStat(t, simKern1.decay, simKern2.decay, ...
                         simKern1.delay, simKern2.delay, sigma);
    h2 = simXsimComputeDiagHStat(t, simKern2.decay, simKern1.decay, ...
                         simKern2.delay, simKern1.delay, sigma);
end

sK = 0.5 * (h1 + h2);
if (isSim1Normalised == false) || (isSim2Normalised == false)
  sK = sqrt(pi) * sigma * sK;
end

if isfield(simKern1, 'isVarS') && (simKern1.isVarS)
    k = sK;
else
    if isfield(simKern1, 'isNegativeS') && (simKern1.isNegativeS == true)
        k = simKern1.sensitivity * sK;
    else
        k = sqrt(simKern1.variance) *sK;
    end
    if isfield(simKern2, 'isNegativeS') && (simKern2.isNegativeS == true)
        k = simKern2.sensitivity * k;
    else
        k = sqrt(simKern2.variance) * k;
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function tstr=mktstr
%	return time string in hh:mm:ss format
%	by Hideki Kawahara
%	05/Jan./1995

anatime=fix(clock);
tstr=[num2str(anatime(4)) ':' num2str(anatime(5)) ':' num2str(anatime(6))];

author: HidekiKawahara


text:
function dG0_prime = Transform(pseudoisomers, pH, I, T)
% Calculate pseudoisomer group standard transformed Gibbs energy of
% formation at specified pH, ionic strength and temperature.
%
% USAGE:
%
%    dG0_prime = Transform(pseudoisomers, pH, I, T)
%
% INPUTS:
%    pseudoisomers:    `p x 3` matrix with a row for each of the p pseudoisomers
%                      in the group, and the following columns:
%
%                        1. Standard Gibbs energy of formation,
%                        2. Number of hydrogen atoms,
%                        3. Charge.
%    pH:               pH.
%    I:                Ionic strength in mol/L.
%    T:                Temperature in Kelvin.
%
% OUTPUT:
%    dG0_prime:        Pseudoisomer group standard transformed Gibbs energy of
%                      formation in kJ/mol.
%
% .. Authors:
%       - Elad Noor, Nov. 2012
%       - Hulda SH, Nov. 2012, Added temperature dependent alpha.

R = 8.3144621e-3; % Gas constant in kJ/(K*mol)
alpha = (9.20483*T)/10^3 - (1.284668*T^2)/10^5 + (4.95199*T^3)/10^8; % Approximation of the temperature dependency of ionic strength effects
DH = (alpha * sqrt(I)) / (1 + 1.6 * sqrt(I)); % Debye Huckel

% dG0' = dG0 + nH * (RTlog(10) * pH + DH) + charge^2 * DH;
dG0_prime_vector = pseudoisomers(:, 1) + ...
    pseudoisomers(:, 2) * (R*T*log(10)*pH + DH) - ...
    pseudoisomers(:, 3).^2 * DH;

try
    dG0_prime = -R * T * maxstar(dG0_prime_vector / (-R * T));
catch
    disp(dG0_prime_vector)
    return
end

author: opencobra

repo: cobratoolbox


text:
function sparse_grid_gl_test01 ( dim_min, dim_max, level_max_min, ...
  level_max_max )

%*****************************************************************************80
%
%% TEST01 tests SPARSE_GRID_GL_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_MIN, the minimum spatial dimension to consider.
%
%    Input, integer DIM_MAX, the maximum spatial dimension to consider.
%
%    Input, integer LEVEL_MAX_MIN, the minimum value of LEVEL_MAX to consider.
%
%    Input, integer LEVEL_MAX_MAX, the maximum value of LEVEL_MAX to consider.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  SPARSE_GRID_GL_SIZE returns the number of distinct\n' );
  fprintf ( 1, '  points in a Gauss-Legendre sparse grid.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Note that, unlike most sparse grids, a sparse grid\n' );
  fprintf ( 1, '  based on Gauss-Legendre points is almost entirely \n' );
  fprintf ( 1, '  NOT nested.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Hence the point counts should be much higher than for a grid of\n' );
  fprintf ( 1, '  the same level, but using rules such as Fejer1 or Fejer2 or\n' );
  fprintf ( 1, '  Gauss-Patterson or Newton-Cotes-Open or Newton-Cotes-Open-Half.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Each sparse grid is of spatial dimension DIM,\n' );
  fprintf ( 1, '  and is made up of all product grids of levels up to LEVEL_MAX.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   DIM: ' );

  for dim_num = dim_min : dim_max
    fprintf ( 1, '  %8d', dim_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   LEVEL_MAX\n' );
  fprintf ( 1, '\n' );

  for level_max = level_max_min : level_max_max
    fprintf ( 1, '    %4d', level_max );
    for dim_num = dim_min : dim_max
      point_num = sparse_grid_gl_size ( dim_num, level_max );
      fprintf ( 1, '  %8d', point_num );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rez = fitTemplates(rez, DATA, uproj)

nt0             = rez.ops.nt0;
rez.ops.nt0min  = ceil(20 * nt0/61);

ops = rez.ops;

rng('default');
rng(1);

Nbatch      = rez.temp.Nbatch;
Nbatch_buff = rez.temp.Nbatch_buff;

Nfilt 	= ops.Nfilt; %256+128;

ntbuff  = ops.ntbuff;
NT  	= ops.NT;

Nrank   = ops.Nrank;
Th 		= ops.Th;
maxFR 	= ops.maxFR;

Nchan 	= ops.Nchan;

batchstart = 0:NT:NT*(Nbatch-Nbatch_buff);

delta = NaN * ones(Nbatch, 1);
iperm = randperm(Nbatch);

switch ops.initialize
    case 'fromData'
        WUinit = optimizePeaks(ops,uproj);%does a scaled kmeans
        dWU    = WUinit(:,:,1:Nfilt);
        %             dWU = alignWU(dWU);
    otherwise
        if ~isempty(getOr(ops, 'initFilePath', [])) && ~getOr(ops, 'saveInitTemps', 0)
            load(ops.initFilePath);
            dWU = WUinit(:,:,1:Nfilt);
        else
            initialize_waves0;

            ipck = randperm(size(Winit,2), Nfilt);
            W = [];
            U = [];
            for i = 1:Nrank
                W = cat(3, W, Winit(:, ipck)/Nrank);
                U = cat(3, U, Uinit(:, ipck));
            end
            W = alignW(W, ops);

            dWU = zeros(nt0, Nchan, Nfilt, 'single');
            for k = 1:Nfilt
                wu = squeeze(W(:,k,:)) * squeeze(U(:,k,:))';
                newnorm = sum(wu(:).^2).^.5;
                W(:,k,:) = W(:,k,:)/newnorm;

                dWU(:,:,k) = 10 * wu;
            end
            WUinit = dWU;
        end
end
if getOr(ops, 'saveInitTemps', 0)
    if ~isempty(getOr(ops, 'initFilePath', []))
        save(ops.initFilePath, 'WUinit')
    else
       warning('cannot save initialization templates because a savepath was not specified in ops.saveInitTemps');
    end
end

[W, U, mu, UtU, nu] = decompose_dWU(ops, dWU, Nrank, rez.ops.kcoords);
W0 = W;
W0(NT, 1) = 0;
fW = fft(W0, [], 1);
fW = conj(fW);

nspikes = zeros(Nfilt, Nbatch);
lam =  ones(Nfilt, 1, 'single');

freqUpdate = 100 * 4;
iUpdate = 1:freqUpdate:Nbatch;

dbins = zeros(100, Nfilt);
dsum = 0;
miniorder = repmat(iperm, 1, ops.nfullpasses);
%     miniorder = repmat([1:Nbatch Nbatch:-1:1], 1, ops.nfullpasses/2);

i = 1; % first iteration

epu = ops.epu;

%%
% pmi = exp(-1./exp(linspace(log(ops.momentum(1)), log(ops.momentum(2)), Nbatch*ops.nannealpasses)));
pmi = exp(-1./linspace(1/ops.momentum(1), 1/ops.momentum(2), Nbatch*ops.nannealpasses));
% pmi = exp(-linspace(ops.momentum(1), ops.momentum(2), Nbatch*ops.nannealpasses));

% pmi  = linspace(ops.momentum(1), ops.momentum(2), Nbatch*ops.nannealpasses);
Thi  = linspace(ops.Th(1),                 ops.Th(2), Nbatch*ops.nannealpasses);
if ops.lam(1)==0
    lami = linspace(ops.lam(1), ops.lam(2), Nbatch*ops.nannealpasses);
else
    lami = exp(linspace(log(ops.lam(1)), log(ops.lam(2)), Nbatch*ops.nannealpasses));
end

if Nbatch_buff<Nbatch
    fid = fopen(ops.fproc, 'r');
end

nswitch = [0];
msg = [];
fprintf('Time %3.0fs. Optimizing templates ...\n', toc)
while (i<=Nbatch * ops.nfullpasses+1)
    % set the annealing parameters
    if i<Nbatch*ops.nannealpasses
        Th      = Thi(i);
        lam(:)  = lami(i);
        pm      = pmi(i);
    end

    % some of the parameters change with iteration number
    Params = double([NT Nfilt Th maxFR 10 Nchan Nrank pm epu nt0]);

    % update the parameters every freqUpdate iterations
    if i>1 &&  ismember(rem(i,Nbatch), iUpdate) %&& i>Nbatch
        dWU = gather_try(dWU);

        % break bimodal clusters and remove low variance clusters
        if  ops.shuffle_clusters &&...
                i>Nbatch && rem(rem(i,Nbatch), 4*400)==1    % i<Nbatch*ops.nannealpasses
            [dWU, dbins, nswitch, nspikes, iswitch] = ...
                replace_clusters(dWU, dbins,  Nbatch, ops.mergeT, ops.splitT, WUinit, nspikes);
        end

        dWU = alignWU(dWU, ops);

        % restrict spikes to their peak group
        %         dWU = decompose_dWU(dWU, kcoords);

        % parameter update
        [W, U, mu, UtU, nu] = decompose_dWU(ops, dWU, Nrank, rez.ops.kcoords);

        if ops.GPU
            dWU = gpuArray(dWU);
        else
            W0 = W;
            W0(NT, 1) = 0;
            fW = fft(W0, [], 1);
            fW = conj(fW);
        end

        NSP = sum(nspikes,2);
        if ops.showfigures
%             set(0,'DefaultFigureWindowStyle','docked')
%             figure;
            subplot(2,2,1)
            for j = 1:10:Nfilt
                if j+9>Nfilt;
                    j = Nfilt -9;
                end
                plot(log(1+NSP(j + [0:1:9])), mu(j+ [0:1:9]), 'o');
                xlabel('log of number of spikes')
                ylabel('amplitude of template')
                hold all
            end
            axis tight;
            title(sprintf('%d  ', nswitch));
            subplot(2,2,2)
            plot(W(:,:,1))
            title('timecourses of top PC')

            subplot(2,2,3)
            imagesc(U(:,:,1))
            title('spatial mask of top PC')

            drawnow
        end
        % break if last iteration reached
        if i>Nbatch * ops.nfullpasses; break; end

        % record the error function for this iteration
        rez.errall(ceil(i/freqUpdate))          = nanmean(delta);

    end

    % select batch and load from RAM or disk
    ibatch = miniorder(i);
    if ibatch>Nbatch_buff
        offset = 2 * ops.Nchan*batchstart(ibatch-Nbatch_buff);
        fseek(fid, offset, 'bof');
        dat = fread(fid, [NT ops.Nchan], '*int16');
    else
        dat = DATA(:,:,ibatch);
    end

    % move data to GPU and scale it
    if ops.GPU
        dataRAW = gpuArray(dat);
    else
        dataRAW = dat;
    end
    dataRAW = single(dataRAW);
    dataRAW = dataRAW / ops.scaleproc;

    % project data in low-dim space
    data = dataRAW * U(:,:);

    if ops.GPU
        % run GPU code to get spike times and coefficients
        [dWU, ~, id, x,Cost, nsp] = ...
            mexMPregMU(Params,dataRAW,W,data,UtU,mu, lam .* (20./mu).^2, dWU, nu);
    else
        [dWU, ~, id, x,Cost, nsp] = ...
            mexMPregMUcpu(Params,dataRAW,fW,data,UtU,mu, lam .* (20./mu).^2, dWU, nu, ops);
    end

    dbins = .9975 * dbins;  % this is a hard-coded forgetting factor, needs to become an option
    if ~isempty(id)
        % compute numbers of spikes
        nsp                = gather_try(nsp(:));
        nspikes(:, ibatch) = nsp;

        % bin the amplitudes of the spikes
        xround = min(max(1, int32(x)), 100);

        dbins(xround + id * size(dbins,1)) = dbins(xround + id * size(dbins,1)) + 1;

        % estimate cost function at this time step
        delta(ibatch) = sum(Cost)/1e3;
    end

    % update status
    if ops.verbose  && rem(i,20)==1
        nsort = sort(round(sum(nspikes,2)), 'descend');
        fprintf(repmat('\b', 1, numel(msg)));
        msg = sprintf('Time %2.2f, batch %d/%d, mu %2.2f, neg-err %2.6f, NTOT %d, n100 %d, n200 %d, n300 %d, n400 %d\n', ...
            toc, i,Nbatch* ops.nfullpasses,nanmean(mu(:)), nanmean(delta), round(sum(nsort)), ...
            nsort(min(size(W,2), 100)), nsort(min(size(W,2), 200)), ...
            nsort(min(size(W,2), 300)), nsort(min(size(W,2), 400)));
        fprintf(msg);
    end

    % increase iteration counter
    i = i+1;
end

% close the data file if it has been used
if Nbatch_buff<Nbatch
    fclose(fid);
end

if ~ops.GPU
   rez.fW = fW; % save fourier space templates if on CPU
end
rez.dWU               = gather_try(dWU);
rez.nspikes               = nspikes;
% %%

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
function S = globMatrixIFE3DPreCond(fun,coef,mesh,femI,fem1,fem2,dindx,dindy)

%% USAGE: generate stiffness global matrix on a 3D mesh
%
% INPUTS:
% fun --- coefficient function
% mesh --- a struct data contains very rich mesh information.
% fem1 --- global DoF for test function space
% fem2 --- global DoF for trial function space
%
% OUTPUTS:
% [IN JN XN] --- triplets of the sparse matrix from regular elements.
% [II JI XI] --- triplets of the sparse matrix from interface elements.

% Last Modified: 08/07/2020 by Xu Zhang

%% 0. Initializaiton
if strcmp(fem1.type,'P1')||strcmp(fem1.type,'DGP1')||strcmp(fem1.type,'CR')
    feEvalBas1 = @evalP1Bas3D;
elseif strcmp(fem1.type,'P2')||strcmp(fem1.type,'DGP2')
    feEvalBas1 = @evalP2Bas3D;
end

if strcmp(fem2.type,'P1')||strcmp(fem2.type,'DGP1')||strcmp(fem2.type,'CR')
    feEvalBas2 = @evalP1Bas3D;
elseif strcmp(fem2.type,'P2')||strcmp(fem2.type,'DGP2')
    feEvalBas2 = @evalP2Bas3D;
end

%% 1. Matrix on noninterface elements
dof1 = fem1.ldof; dof2 = fem2.ldof; nloc = dof1*dof2;
ntID = find(mesh.tLoc > 0); ntN = length(ntID);
AN = fem1.area(ntID);
gxN = fem1.gx(ntID,:); gyN = fem1.gy(ntID,:); gzN = fem1.gz(ntID,:); gw = fem1.gw;
XN = zeros(nloc*ntN, 1);

coefN = feval(fun,gxN,gyN,gzN);
Ibas = cell(dof1,1);
Jbas = cell(dof2,1);
for i = 1:dof1
    Ibas{i} = feEvalBas1(fem1.bas(ntID,:,i), gxN, gyN, gzN, dindx);
end
for j = 1:dof2
    Jbas{j} = feEvalBas2(fem2.bas(ntID,:,j), gxN, gyN, gzN, dindy);
end

IN = reshape(repmat(fem1.t(ntID,:),4,1),nloc*ntN,1);
JN = repmat(reshape(fem2.t(ntID,:),dof2*ntN,1),4,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XN(ind+1:ind+ntN) = AN.*(sum(((Ibas{i}.*(coefN.*Jbas{j})).*gw'),2));
        ind = ind + ntN;
    end
end
ID = find(XN~=0);
SN = sparse(IN(ID),JN(ID),XN(ID),size(fem1.p,1),size(fem2.p,1));

%% 2. Matrix on interface elements
AI = femI.area; gw = femI.gw; gxI = femI.gx; gyI = femI.gy; gzI = femI.gz;
ntI = size(femI.t,1); % not number of interface element, but quadrature element
XI = zeros(nloc*ntI, 1);

coefI = (coef(1)+coef(2))/2;
Ibas = cell(dof1,1);
Jbas = cell(dof2,1);
for i = 1:dof1
    Ibas{i} = feEvalBas1(femI.bas(:,:,i), gxI, gyI, gzI, dindx);
end
for j = 1:dof2
    Jbas{j} = feEvalBas2(femI.bas(:,:,j), gxI, gyI, gzI, dindy);
end

II = reshape(repmat(femI.t,4,1),nloc*ntI,1);
JI = repmat(reshape(femI.t,dof2*ntI,1),4,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XI(ind+1:ind+ntI) = AI.*(sum(((Ibas{i}.*(coefI.*Jbas{j})).*gw'),2));
        ind = ind + ntI;
    end
end
SI = sparse(II,JI,XI,size(fem1.p,1),size(fem2.p,1));
S = SN + SI;

repo: ifem


text:
%#  NLP written by GAMS Convert at 06/20/02 11:29:54
%#
%#  Equation counts
%#     Total       E       G       L       N       X
%#         3       2       1       0       0       0
%#
%#  Variable counts
%#                 x       b       i     s1s     s2s      sc      si
%#     Total    cont  binary integer    sos1    sos2   scont    sint
%#         5       5       0       0       0       0       0       0
%#  FX     0       0       0       0       0       0       0       0
%#
%#  Nonzero counts
%#     Total   const      NL     DLL
%#        12       6       6       0
%#
%#  Reformualtion has removed 1 variable and 1 equation
%
%
%var x1 := 50, >= 50, <= 200;
%var x2 := 37.5, >= 37.5, <= 150;
%var x3 := 45, >= 45, <= 180;
%var x4;
%
%minimize obj: 0.00533*x1^2 + 11.669*x1 + 0.00889*x2^2 + 10.333*x2 + 0.00741*x3^
%              2 + 10.833*x3 + 653.1;
%
%subject to
%
%e2:  - (0.01*(0.0676*x1*x1 + 0.00953*x1*x2 - 0.00507*x1*x3 + 0.00953*x2*x1 +
%    0.0521*x2*x2 + 0.00901*x2*x3 - 0.00507*x3*x1 + 0.00901*x3*x2 + 0.0294*x3*x3
%    ) - 0.000766*x1 - 3.42e-5*x2 + 0.000189*x3) + x4 = 0.040357;
%
%e3:    x1 + x2 + x3 - x4 >= 210;

function test_ipopt

  auxdata = {} ;

  options.lb = [ 50, 37.5, 45, -Inf ] ;  % Lower bound on the variables.
  options.ub = [ 200, 150, 180, Inf ] ;  % Upper bound on the variables.

  % The constraint functions are bounded to zero
  options.cl = [ 0, 0 ]; %  constraints
  options.cu = [ 0, Inf ];

  % Set up the auxiliary data.
  options.auxdata = auxdata ;

  % Set the IPOPT options.
  options.ipopt.jac_d_constant   = 'no';
  options.ipopt.hessian_constant = 'no';
  options.ipopt.mu_strategy      = 'adaptive';
  options.ipopt.max_iter         = 400;
  options.ipopt.tol              = 1e-10;

  % The callback functions.
  funcs.objective         = @objective;
  funcs.constraints       = @constraints;
  funcs.gradient          = @gradient;
  funcs.jacobian          = @jacobian;
  funcs.jacobianstructure = @jacobianstructure;
  if true
    funcs.hessian           = @hessian;
    funcs.hessianstructure  = @hessianstructure;
    options.ipopt.derivative_test = 'second-order';
  else
    options.ipopt.hessian_approximation      = 'limited-memory';
    %options.ipopt.limited_memory_update_type = 'bfgs' ; % {bfgs}, sr1 = 6; % {6}
    %options.ipopt.limited_memory_update_type = 'sr1' ;
    options.ipopt.limited_memory_update_type = 'bfgs' ; % {bfgs}, sr1 = 6; % {6}
  end

  % Run IPOPT.
  x0 = [50, 37.5, 45, 0] ;

  tic
  [x, info] = ipopt_auxdata(x0,funcs,options);
  elapsed = toc ;

  info;

  x

end

%%
% map the indices with the corresponding index in the spase matrix
function f = objective(x,auxdata)
  f = 0.00533*x(1)^2 + 11.669*x(1) + ...
      0.00889*x(2)^2 + 10.333*x(2) + ...
      0.00741*x(3)^2 + 10.833*x(3) + 653.1 ;
end

%%
% map the indices with the corresponding index in the spase matrix
function g = gradient(x,auxdata)
  g = [ 0.01066*x(1) + 11.669, 0.01778*x(2) + 10.333, 0.01482*x(3) + 10.833, 0 ] ;
end

function f = constraints(x,auxdata)
  f = zeros(2,1) ;
  f(1) = - (0.01*(0.0676*x(1)^2 + 0.00953*x(1)*x(2) - 0.00507*x(1)*x(3) + ...
                  0.00953*x(2)*x(1) + 0.0521*x(2)^2 + 0.00901*x(2)*x(3) - ...
                  0.00507*x(3)*x(1) + 0.00901*x(3)*x(2) + 0.0294*x(3)*x(3) ) ...
         - 0.000766*x(1) - 3.42e-5*x(2) + 0.000189*x(3)) + x(4) - 0.040357 ; % = 0
  f(2) = x(1) + x(2) + x(3) - x(4) - 210 ; % >= 0
end

function jac = jacobian(x,auxdata)
  jac = [ -0.001352*x(1) - 0.0001906*x(2) + 0.0001014*x(3) + 0.000766, ...
          -0.0001906*x(1) - 0.001042*x(2) - 0.0001802*x(3) + 0.0000342, ...
           0.0001014*x(1) - 0.0001802*x(2) - 0.000588*x(3) - 0.000189, ...
           1 ; 1, 1, 1, -1 ] ;
  jac = sparse(jac) ;
end

function jac = jacobianstructure(auxdata)
  jac = sparse(ones(2,4)) ;
end

function H = hessian(x, sigma, lambda, auxdata)
  H1 = [ 0.01066 0       0       0 ; ...
         0       0.01778 0       0 ; ...
         0       0       0.01482 0 ; ...
         0       0       0       0 ] ;
  H2 = [ -0.1352e-2,          0,         0, 0 ; ...
         -0.1906e-3, -0.1042e-2,         0, 0 ; ...
          0.0001014, -0.1802e-3, -0.588e-3, 0 ; ...
          0,                  0,         0, 0 ] ;
  H = sparse(sigma*H1 + lambda(1)*H2) ;
end

function H = hessianstructure(auxdata)
  H = sparse([ 1 0 0 0 ; 1 1 0 0 ; 1 1 1 0 ; 1 1 1 1 ]) ;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function y = cvx_s_symmetric( m, n )
%CVX_S_SYMMETRIC Symmetric matrices (lower triangle storage).

if m ~= n,
    error( 'Symmetric structure requires square matrices.' );
end

nsq = n * n;
ntr = 0.5 * ( nsq + n );
c  = 0 : n - 1;
c  = c( ones( 1, n ), : );
r  = c';
mn = min( r, c );
mx = max( r, c );
y  = mx + mn .* ( n - 0.5 * ( mn + 1 ) ) + 1;
y  = sparse( y( : ), 1 : nsq, 1, ntr, nsq );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
clear all; close all;
I=imread('glass.png');
J=imcomplement(I);
figure;
subplot(121);
imshow(uint8(I));
subplot(122);
imshow(uint8(J));

author: HuangCongQing


text:
%% vecPair2Rot
% Below is a demonstration of the features of the |vecPair2Rot| function

%%
clear; close all; clc;

%% Syntax
% |[R]=vecPair2Rot(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.


text:
function extended_rosenbrock_test ( )

%*****************************************************************************80
%
%% EXTENDED_ROSENBROCK_TEST works with the extended Rosenbrock function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXTENDED_ROSENBROCK_TEST:\n' );
  fprintf ( 1, '  Test COMPASS_SEARCH with the extended Rosenbrock function.\n' );
  m = 4;
  delta_tol = 0.00001;
  delta = 0.3;
  k_max = 20000;

  x = [ - 1.2, 1.0,  -1.5, 1.2 ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', extended_rosenbrock ( m, x ) );

  [ x, fx, k ] = compass_search ( @extended_rosenbrock, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );
%
%  Demonstrate correct minimizer.
%
  x = [ 1.0, 1.0, 1.0, 1.0 ];
  r8vec_print ( m, x, '  Correct minimizer X*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', extended_rosenbrock ( m, x ) );

  return
end


text:
function [theCombo,recurData]=getNextCapacitatedCombo(param1,numBins,totalItems)
%%GETNEXTCAPACITATEDCOMBO Get the next combination of placing totalItems
%               items into numBins unlabeled bins where each bin can hold
%               numPerBin items. If numPerBin=1, then this is just
%               generating standard combinations of numBins items taken
%               from a total of totalItems. The capacitated combinations
%               are equivalent to putting totalItems labeled balls into
%               numBins unlabeled bins such that each bin holds numPerBin
%               items and the ordering of the items in the bins does not
%               matter.
%
%INPUTS: param1 If the first combination is desired, then more than one
%               input must be given and param1=numPerBin, the number of
%               items in each bin. Otherwise, param1=recurData, the second
%               output of the function from the last call.
%       numBins This input should only be provided if the first capacitated
%               combination in the sequence is desired. This is the number
%               of bins (of capacity numPerBin) into which things can be
%               placed.
%    totalItems This input can only be given when the first combination is
%               desired. This is the total number of items that can be
%               assigned. If numBins is provided, but this is omitted or an
%               empty matrix is passed, then a value of numBins*numPerBin
%               is used, the minimum allowed value.
%
%OUTPUTS: theCombo The numPerBinXnumBins next capacitated combination. If
%                  the final capacitated combination has been passed, then
%                  an empty matrix is returned.
%        recurData A data structure that can be passed back to this
%                  function to get the next combination.
%
%The algorithm is the same as that used in genAllCapacitatedCombos. See the
%comments to genAllCapacitatedCombos for more information. The total number
%of capaciatted combinations is given by numCapacitatedCombos.
%
%EXAMPLE:
%This example just shows that the capacitated combinations produced by this
%function come in the same order as those produced by the function
%genAllCapacitatedCombos. When run, it collects the combinations by
%repeatedly calling this function and then asserts that they are all equal
%what genAllCapacitatedCombos returns. The result shoud be true, so there
%should be no error.
% numPerBin=2;
% numBins=4;
% totalItems=4*2+2;
%
% totalCapCombos=numCapacitatedCombos(numPerBin,numBins,totalItems);
% allCapCombos=zeros(numPerBin,numBins,totalCapCombos);
%
% [theCombo,recurData]=getNextCapacitatedCombo(numPerBin,numBins,totalItems);
% allCapCombos(:,:,1)=theCombo;
%
% for comboIdx=2:totalCapCombos
%     [theCombo,recurData]=getNextCapacitatedCombo(recurData);
%     allCapCombos(:,:,comboIdx)=theCombo;
% end
%
% allBatchCombos=genAllCapacitatedCombos(numPerBin,numBins,totalItems);
% assert(all(allCapCombos(:)==allBatchCombos(:)))
%The assertion should not raise an error.
%
%July 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin>=2)
    if(nargin<3||isempty(totalItems))
        totalItems=param1*numBins;
    end

    %If this is the first capacitated combination, then create the
    %recurData structure for the first time.
    numPerBin=param1;

    %The total number of items to assign in combinations across the bins.
    n=numBins*numPerBin;

    %Save the values in recurData that do not change with the iterations.
    recurData=[];
    recurData.numPerBin=numPerBin;
    recurData.numBins=numBins;
    recurData.n=n;
    recurData.totalItems=totalItems;

    %Special case for standard combinations.
    if(numPerBin==1)
        curIdxCombos=1:numBins;
        recurData.curIdxCombos=curIdxCombos;
        theCombo=curIdxCombos;
        recurData.curBin=1;
        recurData.curSet=1;
        recurData.curCombos=[];
        recurData.idxLists=[];
        recurData.backtracking=[];

        I=1:n;
        recurData.I=getNextCombo(I,totalItems,1);
        return
    end

    curIdxCombos=zeros(numPerBin,numBins);
    curCombos=zeros(numPerBin-1,numBins);

    %This holds the indices that can be assigned in each bin and in subsequent
    %bins.
    idxLists=zeros(n,numBins);
    idxLists(:,1)=1:n;

    curSet=1;
    curBin=1;
    backtracking=false;
    %When n~=totalItems, after each base capacitated combination is
    %generated, we have to go through all combinations of indices that make
    %up the combinations. This is the current index combination. When
    %empty, it means that we haven't gotten a base capacitated combination
    %yet.
    I=[];
    recurData.I=I;
else
    recurData=param1;

    if(isempty(recurData))
       theCombo=[];
       return;
    end

    %Extract the values that change on each iteration.
    curIdxCombos=recurData.curIdxCombos;
    curCombos=recurData.curCombos;
    idxLists=recurData.idxLists;
    curSet=recurData.curSet;
    curBin=recurData.curBin;
    backtracking=recurData.backtracking;
    I=recurData.I;
end

%Get the values that do not change between calls.
numPerBin=recurData.numPerBin;
numBins=recurData.numBins;
n=recurData.n;
totalItems=recurData.totalItems;
numPerBinM1=numPerBin-1;

if(~isempty(I))
    %Just get the next combination of indices to the capcitated combination
    %and return that.
    theCombo=reshape(I(curIdxCombos),[numPerBin,numBins]);

    %If past the final combination, then I is an empty matrix and it will
    %know to go onto the next capacitated base combination.
    I=getNextCombo(I,totalItems,1);
    recurData.I=I;
    return
elseif(numPerBin==1)
    %The special case of standard combos.
    theCombo=[];
    return;
end

while(curBin>0)
    %This is the number of elements that can be assigned in the current
    %bin including making a hard assignment on the first index.
    numInIdx=n-(curBin-1)*numPerBin;

    if(backtracking==false)
        %We are just entering this bin.
        %The first element in each bin is fixed to the first index in the
        %idxList.
        curIdxCombos(1,curBin)=idxLists(1,curBin);

        %In this level, we have to go through all combinations of the
        %remaining n-curBin*(numPerBin-1); unassigned elements of
        %idxLists(:,curBin) in this bin. Each time an assignment is made in
        %this bin, those elements are removed from idxLists going into the
        %next higher level so that we never get repeats.
        curIdxCombos(2:numPerBin,curBin)=idxLists(2:numPerBin,curBin);
        %Current combination of the indices in idxLists in this bin.
        curCombos(:,curBin)=1:numPerBinM1;

        if(curBin==numBins)
            %If we are in the final bin, then save the combo and move to
            %the previous bin. In practice, subsequent calls to this
            %function will nto go to the previous bin until all
            %combinations of index assignments (via I) are assigned.
            I=1:n;
            theCombo=curIdxCombos;
            curSet=curSet+1;

            backtracking=true;
            curBin=curBin-1;

            %Save the values that change during iteration to recurData
            %before returning.
            recurData.curIdxCombos=curIdxCombos;
            recurData.curCombos=curCombos;
            recurData.idxLists=idxLists;
            recurData.curSet=curSet;
            recurData.curBin=curBin;
            recurData.backtracking=backtracking;
            recurData.I=getNextCombo(I,totalItems,1);
            return;
        else
            %Otherwise, set idxLists for the next level and move to the
            %next bin.numInIdxNext
            numInIdxNext=numInIdx-numPerBin;
            idxLists(1:numInIdxNext,curBin+1)=idxLists((numPerBin+1):numInIdx,curBin);

            curBin=curBin+1;
        end
    else
        %If here, we are backtracking.

        %We subtract 1 because we have made a hard assignment on the first
        %index.
        nextCombo=getNextCombo(curCombos(:,curBin),numInIdx-1,1);
        if(isempty(nextCombo))
            %Continue backtracking.
            curBin=curBin-1;
            continue;
        end
        %Otherwise, record the current combination.
        curCombos(:,curBin)=nextCombo;

        %Record the indices for the current combination.
        nextCombo=nextCombo+1;
        curIdxCombos(2:numPerBin,curBin)=idxLists(nextCombo,curBin);

        %Now, set the values in idxLists at the next bin level to those
        %that are in idxLists but not in nextCombo. We note that the
        %indices are given in increasing order and that the first index is
        %skipped because we know it is assigned.
        curAssignedIdx=2;
        curNextIdx=1;
        for curListIdx=2:numInIdx
            if(idxLists(curListIdx,curBin)==curIdxCombos(curAssignedIdx,curBin))
                curAssignedIdx=curAssignedIdx+1;
                if(curAssignedIdx>numPerBin)
                    %Then all of the rest of the indices are assigned.
                    numInIdxNext=numInIdx-numPerBin;
                    idxLists(curNextIdx:numInIdxNext,curBin+1)=idxLists((curListIdx+1):numInIdx,curBin);
                    break;
                end
            else
                idxLists(curNextIdx,curBin+1)=idxLists(curListIdx,curBin);
                curNextIdx=curNextIdx+1;
            end
        end

        %Go to the next bin.
        curBin=curBin+1;
        backtracking=false;
    end
end

%If we get here, we passed the final combiantion.
theCombo=[];
recurData=[];
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:
% C = FindClustermeans(gIX,M);
% [C_trialAvr,C_trialRes,C_score,C_d2var_perstim] = GetTrialAvrLongTrace(hGUI,C);

%% Init load
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%% Load Fish8
i_fish = 8;
[cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,[1,1]);

fishset = getappdata(hfig,'fishset');
stim = getappdata(hfig,'stim');
behavior = getappdata(hfig,'behavior');
fpsec = getappdata(hfig,'fpsec');

%% crop
cropIX = 3001:3600;
stim_crop = stim(cropIX);
behavior_crop = behavior(:,cropIX);%(2,3001:3600);

% stim z
regressors_s = GetStimRegressor(stim_crop,fishset,i_fish);
reg_s = regressors_s(8).im;
% motor
regressors_m = GetMotorRegressor(behavior_crop);
reg_m = regressors_m(2).im; %1
reg_motor_n = reg_m/sqrt(sum(reg_m.*reg_m));
%%
reg_stim_n = reg_s/sqrt(sum(reg_s.*reg_s));
%%

%% get tRes/tAvr
C = reg_motor_n;
% GetTrialAvrLongTrace
period = 150; % 120 for PT, 150 for OMR
C_3D_0 = reshape(C,size(C,1),period,[]);

%     C_3D = zscore(C_3D_0,0,2);
%     C_d2var_perstim = nanmean(nanstd(C_3D,0,3),2);
%     C_score = C_d2var_perstim;

C_period = mean(C_3D_0,3);%prctile(C_3D_0,20,3);%mean(C_3D_0,3);
nPeriods = round(size(C,2)/period);
C_trialAvr = repmat(C_period,1,nPeriods);

C_trialRes = C-C_trialAvr;

%% tRes/tAvr demo trace plot
xv = (1:600)/fpsec;

% find vertical stim time-references
x_stim = stim_crop';
x_stim(x_stim~=11) = 9;
xlines = xv(find(diff([0;x_stim])~=0));

figure('Position',[50,200,300,200]); hold on

% plot vertical reference shadings
for i = 2:2:length(xlines)
    x_ = [xlines(i),xlines(i+1)];
    y_ = [-1,0.2];
    x = [x_(1),x_(2),x_(2),x_(1)];
    y = [y_(1),y_(1),y_(2),y_(2)];
   patch(x, y, [1,0.9,0.9],'Edgecolor','w');
end

% plot traces
plot(xv,reg_stim_n,'k')

yoffset = -0.3;
plot(xv,reg_motor_n+yoffset,'k')

yoffset = -0.6;
plot(xv,C_trialAvr+yoffset,'k')

yoffset = -0.9;
plot(xv,C_trialRes+yoffset,'k')

xlim([0,300])
% xlabel('Time (sec)')
ax = gca;
ax.YTick = [-0.9,-0.6,-0.3,0];
ax.YTickLabel = {'trial-Res.','trial-Avr.','motor','stim.'};

ylim([-1.2,0.2]);
set(gca,'xcolor','w','xtick',[]);
set(gca,'TickLength',[0,0]);
% plot scale bar
plot([xv(1)+5,xv(40)+5],[-1.1,-1.1],'k','linewidth',1.5);
text(xv(1)+5,-1.2,'20 sec')

%% figS3h: L/R demo trace plot
i_fish = 2;
[cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,[1,1]);

fishset = getappdata(hfig,'fishset');
stim = getappdata(hfig,'stim');
behavior = getappdata(hfig,'behavior');
fpsec = getappdata(hfig,'fpsec');

%% crop
cropIX = 3001:3300;%3001:3600;
stim_crop = stim(cropIX);
behavior_crop = behavior(:,cropIX);%(2,3001:3600);
regressors_s = GetStimRegressor(stim_crop,fishset,i_fish);
regressors_m = GetMotorRegressor(behavior_crop);

reg_mL = regressors_m(1).im; %1
reg_mR = regressors_m(2).im; %1
reg_mL_n = reg_mL/sqrt(sum(reg_mL.*reg_mL));
reg_mR_n = reg_mR/sqrt(sum(reg_mR.*reg_mR));
lrAvr = 0.5*(reg_mL_n+reg_mR_n);
lrAvr_n = lrAvr/sqrt(sum(lrAvr.*lrAvr));
lrRes1 = reg_mL_n-lrAvr;
lrRes_n1 = lrRes1/sqrt(sum(lrRes1.*lrRes1));
lrRes2 = reg_mR_n-lrAvr;
lrRes_n2 = lrRes2/sqrt(sum(lrRes2.*lrRes2));
%%
figure('Position',[50,200,400,300]); hold on
xv = (1:length(cropIX))/fpsec;

% plot vertical reference shadings
% for i = 2:2:length(xlines)
%     x_ = [xlines(i),xlines(i+1)];
%     y_ = [-1,0.2];
%     x = [x_(1),x_(2),x_(2),x_(1)];
%     y = [y_(1),y_(1),y_(2),y_(2)];
%    patch(x, y, [1,0.9,0.9],'Edgecolor','w');
% end

% plot traces
% plot(xv,reg_stim_n,'k')
spacing = 0.5;
yoffset = 0;
plot(xv,reg_mL_n+yoffset,'k')

yoffset = yoffset-spacing;
plot(xv,reg_mR_n+yoffset,'k')

yoffset = yoffset-spacing;
plot(xv,lrAvr_n+yoffset,'k')

yoffset = yoffset-spacing;
plot(xv,lrRes_n1+yoffset,'k')

yoffset = yoffset-spacing;
plot(xv,lrRes_n2+yoffset,'k')

% xlim([0,300])
% xlabel('Time (sec)')
ax = gca;
ax.YTick = [-spacing*4,-spacing*3,-spacing*2,-spacing,0];
ax.YTickLabel = {'R-Res.','L-Res.','LR-Avr.','R motor','L motor'};

ylim([-spacing*5,spacing]);
set(gca,'xcolor','w','xtick',[]);
set(gca,'TickLength',[0,0]);
% plot scale bar
plot([xv(1)+5,xv(40)+5],[-spacing*4.6,-spacing*4.6],'k','linewidth',1.5);
text(xv(1)+5,-spacing*4.9,'20 sec')

author: xiuyechen

repo: FishExplorer


text:
function x = pvand ( n, alpha, b )

%*****************************************************************************80
%
%% PVAND solves a Vandermonde system A * x = b.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ake Bjorck, Victor Pereyra,
%    Solution of Vandermonde Systems of Equations,
%    Mathematics of Computation,
%    Volume 24, Number 112, October 1970, pages 893-903.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real ALPHA(N), the parameters that define the matrix.
%    The values should be distinct.
%
%    Input, real B(N), the right hand side of the linear system.
%
%    Output, real X(N), the solution of the linear system.
%
  x(1:n) = b(1:n);

  for k = 1 : n - 1
    for j = n : -1 : k + 1
      x(j) = x(j) - alpha(k) * x(j-1);
    end
  end

  for k = n - 1 : -1 : 1
    for j = k + 1 : n
      x(j) = x(j) / ( alpha(j) - alpha(j-k) );
    end
    for j = k : n - 1
      x(j) = x(j) - x(j+1);
    end
  end

  return
end


text:
function res = trialtag(this, varargin)
% Method for getting/setting trial tag
% FORMAT res = trialtag(this, ind, tag)
%   ind = indices of trials
% The user can put any data here that will be attached to
% the respective trial. This is useful e.g. to make sure the
% relation between regressors and data is not broken when
% removing bad trials or merging files.
% _______________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: trialtag.m 6618 2015-12-01 16:25:38Z spm $

res = getset(this, 'trials', 'tag', varargin{:});

author: spm


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% utl_signmatrix
%
% Goal: create signmatrix
% max_d is max degree
%
% Li Shen
% 01/14/2007 - create

function utl_sgm(max_d)
global sgm;

for d=1:max_d
    n = 2*d+1;
    M = (-1).^(1:n^2);
    sgm{d} = reshape(M,n,n);
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef symm

methods(Static)

    function qmf = quad_mirror_filter(num_vanishing_moments)
        % Generates orthonormal quadrature mirror filter for Symmlets
        switch num_vanishing_moments
            case 4
                qmf = [   -.107148901418  -.041910965125  .703739068656   ...
                1.136658243408  .421234534204   -.140317624179  ...
                -.017824701442  .045570345896                   ];
            case 5
                qmf = [   .038654795955   .041746864422   -.055344186117  ...
                .281990696854   1.023052966894  .896581648380   ...
                .023478923136   -.247951362613  -.029842499869  ...
                .027632152958                                   ];
            case 6
                qmf = [   .021784700327   .004936612372   -.166863215412  ...
                -.068323121587  .694457972958   1.113892783926  ...
                .477904371333   -.102724969862  -.029783751299  ...
                .063250562660   .002499922093   -.011031867509  ];
            case 7
                qmf = [   .003792658534   -.001481225915  -.017870431651  ...
                .043155452582   .096014767936   -.070078291222  ...
                .024665659489   .758162601964   1.085782709814  ...
                .408183939725   -.198056706807  -.152463871896  ...
                .005671342686   .014521394762                   ];
            case 8
                qmf = [   .002672793393   -.000428394300  -.021145686528  ...
                .005386388754   .069490465911   -.038493521263  ...
                -.073462508761  .515398670374   1.099106630537  ...
                .680745347190   -.086653615406  -.202648655286  ...
                .010758611751   .044823623042   -.000766690896  ...
                -.004783458512                                  ];
            case 9
                qmf = [   .001512487309   -.000669141509  -.014515578553  ...
                .012528896242   .087791251554   -.025786445930  ...
                -.270893783503  .049882830959   .873048407349   ...
                1.015259790832  .337658923602   -.077172161097  ...
                .000825140929   .042744433602   -.016303351226  ...
                -.018769396836  .000876502539   .001981193736   ];
            case 10
                qmf = [   .001089170447   .000135245020   -.012220642630  ...
                -.002072363923  .064950924579   .016418869426   ...
                -.225558972234  -.100240215031  .667071338154   ...
                1.088251530500  .542813011213   -.050256540092  ...
                -.045240772218  .070703567550   .008152816799   ...
                -.028786231926  -.001137535314  .006495728375   ...
                .000080661204   -.000649589896                  ];
            otherwise
                error('Unsupported number of vanishing moments');
        end
        % Normalize the coefficients
        qmf = qmf ./ sqrt(2);
    end

    function wave = wavelet_function(j, k, n, num_vanishing_moments)
        % Makes an orthogonal wavelet function
        %
        % Inputs:
        %  j, k - scale and location indices
        %  n - signal length (dyadic)
        qmf = spx.wavelet.symm.quad_mirror_filter(num_vanishing_moments);
        w = zeros(1, n);
        % identify the index of the j-th scale at k-th translation
        index = spx.wavelet.dyad_to_index(j, k);
        w(index) = 1;
        wave = spx.wavelet.transform.inverse_periodized_orthogonal(qmf, w, j);
    end

    function wave = scaling_function(j, k, n, num_vanishing_moments)
        % Makes an orthogonal scaling function
        %
        % Inputs:
        %  j, k - scale and location indices
        %  n - signal length (dyadic)
        qmf = spx.wavelet.symm.quad_mirror_filter(num_vanishing_moments);
        w = zeros(1, n);
        % k-th translate in the coarsest part of the wavelet coefficients
        w(k) = 1;
        wave = spx.wavelet.transform.inverse_periodized_orthogonal(qmf, w, j);
    end

end

end

author: indigits

repo: sparse-plex


text:
function [kernel, S] = dat2Kernel(data, kSize)
% kernel = dat2Kernel(data, kSize,thresh)
%
% Function to perform k-space calibration step for ESPIRiT and create
% k-space kernels. Only works for 2D multi-coil images for now.
%
% Inputs:
%       data - calibration data [kx,ky,coils]
%       kSize - size of kernel (for example kSize=[6,6])
%
% Outputs:
%       kernel - k-space kernels matrix (not cropped), which correspond to
%                the basis vectors of overlapping blocks in k-space
%       S      - (Optional parameter) The singular vectors of the
%                 calibration matrix
%
%
% See also:
%           kernelEig
%
% (c) Michael Lustig 2013

[sx,sy,nc] = size(data);
imSize = [sx,sy] ;

tmp = im2row(data,kSize); [tsx,tsy,tsz] = size(tmp);
A = reshape(tmp,tsx,tsy*tsz);

[U,S,V] = svd(A,'econ');

kernel = reshape(V,kSize(1),kSize(2),nc,size(V,2));
S = diag(S);S = S(:);

author: thomaskuestner


text:
report_this_filefun(mfilename('fullpath'));

% drap a colormap onto topography
l = isnan(tmap);
tmap(l) = -10;

%l = tmap< 0.1;
%tmap(l) = 0;

[lat,lon] = meshgrat(tmap,tmapleg);
%[smap,smapleg] = country2mtx('switzerland',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

% tmap = km2deg(tmap/1);
[X , Y]  = meshgrid(gx,gy);

%sw = interp2(lon0,lat0,smap,lon,lat);
l = re4 < 1.2;
re4(l) = 1.2;
l = re4 > 2.6;
re4(l) = 2.6;

ren = interp2(X,Y,re4,lon,lat);

mi = min(min(ren));
l =  isnan(ren);
ren(l) = mi-20;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[s4_south s3_north],'MapLonLimit',[s2_west s1_east])

l = tmap  == -10 &  ren < mi -10;
tmap(l) = nan;

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',4);
tightmap
view([0 90])
camlight; lighting phong
set(gca,'projection','perspective');

load worldlo
%h = displaym(POline); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',1.7)
h2 = displaym(PPpoint);
h = displaym(PPtext); trimcart(h);

load usahi
[la , lo ] = extractm(stateline,'california');

plotm(la, lo,'w','Linewidth',1.4);
zdatam(handlem('allline'),10000) % keep line on surface
zdatam(handlem('alltext'),10000) % keep line on surface

% pl = plotm(ma(:,2),ma(:,1),'hw');
%  set(pl,'LineWidth',1.5,'MarkerSize',12,...
% 'MarkerFaceColor','w','MarkerEdgeColor','k')
plotm(faults(:,2), faults(:,1),'k','Linewidth',2);
zdatam(handlem('allline'),10000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface

j = hsv;
%j = j(64:-1:1,:);
j = [ [ 0.9 0.9 0.9 ] ; j];
caxis([ 1.15 2.7])

colormap(j);

axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',3);

setm(gca,'mlabellocation',2.5)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',2.5)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','Fontweight','bold','FontSize',16,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.8 0.35 0.01 0.3],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','bold','FontSize',16);
set(gcf,'Inverthardcopy','off');


text:
function r8mat_transpose_test ( )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_TEST tests R8MAT_TRANSPOSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2015
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_TRANSPOSE_TEST\n' );
  fprintf ( 1, '  R8MAT_TRANSPOSE transposes an R8MAT.\n'; );

  a = r8mat_indicator ( m, n );
  r8mat_print ( m, n, a, '  Matrix A:' );

  at = r8mat_transpose ( m, n, a );
  r8mat_print ( n, m, at, '  Transposed matrix At:' );

  return
end


text:
function [Nm] = eV2Nm(eV)
% Convert energy or work from electron volts to newton-meters.
% Chad A. Greene 2012
Nm = eV*1.6021773e-19;

author: Sable

repo: mcbench-benchmarks


text:
function newModel=addMets(model,metsToAdd,copyInfo,prefix)
% addMets
%   Adds metabolites to a model
%
%   model       a model structure
%   metsToAdd   the metabolite structure can have the following fields:
%               mets           cell array with unique strings that
%                              identifies each metabolite (opt, default IDs
%                              of new metabolites are numbered with the
%                              prefix defined below)
%               metNames       cell array with the names of each
%                              metabolite
%               compartments   cell array with the compartment of each
%                              metabolite. Should match model.comps.
%                              If this is a string rather than a cell
%                              array it is assumed that all mets are in
%                              that compartment
%               b              Nx1 or Nx2 matrix with equality constraints
%                              for each metabolite (opt, default 0)
%               unconstrained  vector describing if each metabolite is an
%                              exchange metabolite (1) or not (0) (opt,
%                              default 0)
%               inchis         cell array with InChI strings for each
%                              metabolite (opt, default '')
%               metSmiles      cell array with SMILES strings for each
%                              metabolite (opt, default '')
%               metFormulas    cell array with the formulas for each of
%                              the metabolites (opt, default '')
%               metMiriams     cell array with MIRIAM structures (opt,
%                              default [])
%               metCharges     metabolite charge (opt, default NaN)
%               metNotes       cell array with metabolite notes as strings
%                              (opt, default '')
%   copyInfo    when adding metabolites to a compartment where it previously
%               doesn't exist, the function will copy any available annotation
%               from the metabolite in another compartment (opt, default true)
%   prefix      when metsToAdd.mets is not specified, new metabolite IDs
%               are generated with the prefix specified here. If IDs with
%               the prefix are already used in the model then the
%               numbering will start from the highest existing integer+1
%               (opt, default 'm_')
%
%   newModel     an updated model structure
%
%   NOTE: This function does not make extensive checks about MIRIAM formats,
%   forbidden characters or such.
%
%   Example: If multiple metabolites are added at once, the metMiriams cell
%   array should be defined as (example with ChEBI and KEGG):
%
%   metsToAdd.metMiriams{1} = struct('name',{{'chebi';'kegg.compound'}},...
%       'value',{{'CHEBI:18072';'C11821'}});
%   metsToAdd.metMiriams{2} = struct('name',{{'chebi';'kegg.compound'}},...
%       'value',{{'CHEBI:31132';'C12248'}});
%
%   Usage: newModel=addMets(model,metsToAdd,copyInfo,prefix)

if nargin<3
    copyInfo=true;
end
if nargin<4
    prefix='m_';
else
    prefix=char(prefix);
end

newModel=model;

if isempty(metsToAdd)
    return;
end

%Check some stuff regarding the required fields
if ~isfield(metsToAdd,'mets')
    metsToAdd.mets=generateNewIds(newModel,'mets',prefix,numel(metsToAdd.metNames));
else
    metsToAdd.mets=convertCharArray(metsToAdd.mets);
end
if ~isfield(metsToAdd,'metNames')
    metsToAdd.metNames=metsToAdd.mets;
else
    metsToAdd.metNames=convertCharArray(metsToAdd.metNames);
end
if ~isfield(metsToAdd,'compartments')
    EM='compartments is a required field in metsToAdd';
    dispEM(EM);
else
    metsToAdd.compartments=convertCharArray(metsToAdd.compartments);
    %If only one compartment is given, assume it is for all metabolites
    if numel(metsToAdd.compartments)==1 && numel(metsToAdd.mets)>1
        temp=cell(numel(metsToAdd.mets),1);
        temp(:)=metsToAdd.compartments;
        metsToAdd.compartments=temp;
    end
end

%Number of metabolites
nMets=numel(metsToAdd.mets);
nOldMets=numel(model.mets);
filler=cell(nMets,1);
filler(:)={''};
largeFiller=cell(nOldMets,1);
largeFiller(:)={''};

%Check that no metabolite ids are already present in the model
I=ismember(metsToAdd.mets,model.mets);
if any(I)
    error('One or more elements in metsToAdd.mets are already present in model.mets: %s',...
        strjoin(metsToAdd.mets(I),', '));
else
    newModel.mets=[newModel.mets;metsToAdd.mets(:)];
end

%Check that all the compartments could be found
[I, compMap]=ismember(metsToAdd.compartments,model.comps);
if ~all(I)
    EM='metsToAdd.compartments must match model.comps';
    dispEM(EM);
end

%Check that the metabolite names aren't present in the same compartment.
%Not the neatest way maybe..
t1=strcat(metsToAdd.metNames(:),'***',metsToAdd.compartments(:));
t2=strcat(model.metNames,'***',model.comps(model.metComps));
if any(ismember(t1,t2))
    EM='One or more elements in metsToAdd.metNames already exist in the same compartments as the one it is being added to';
    dispEM(EM);
end

%Some more checks and if they pass then add each field to the structure
if numel(metsToAdd.metNames)~=nMets
    EM='metsToAdd.metNames must have the same number of elements as metsToAdd.mets';
    dispEM(EM);
else
    newModel.metNames=[newModel.metNames;metsToAdd.metNames(:)];
end

if numel(compMap)~=nMets
    EM='metsToAdd.compartments must have the same number of elements as metsToAdd.mets';
    dispEM(EM);
else
    newModel.metComps=[newModel.metComps;compMap(:)];
end

if isfield(metsToAdd,'b')
    if size(metsToAdd.b,1)~=nMets
        EM='metsToAdd.b must have the same number of elements as metsToAdd.mets';
        dispEM(EM);
    else
        %Add empty field if it doesn't exist
        if ~isfield(newModel,'b')
            newModel.b=zeros(nOldMets,1);
        end

        %If the original is only one vector
        if size(metsToAdd.b,2)>size(newModel.b,2)
            newModel.b=[newModel.b newModel.b];
        end
        %Add the new ones
        newModel.b=[newModel.b;metsToAdd.b];
    end
else
    if isfield(newModel,'b')
        %Add the default
        newModel.b=[newModel.b;zeros(nMets,size(newModel.b,2))];
    end
end

if isfield(metsToAdd,'unconstrained')
    if numel(metsToAdd.unconstrained)~=nMets
        EM='metsToAdd.unconstrained must have the same number of elements as metsToAdd.mets';
        dispEM(EM);
    else
        %Add empty field if it doesn't exist
        if ~isfield(newModel,'unconstrained')
            newModel.unconstrained=zeros(nOldMets,1);
        end

        %Add the new ones
        newModel.unconstrained=[newModel.unconstrained;metsToAdd.unconstrained(:)];
    end
else
    if isfield(newModel,'unconstrained')
        %Add the default
        newModel.unconstrained=[newModel.unconstrained;zeros(nMets,1)];
    end
end

if isfield(metsToAdd,'inchis')
    metsToAdd.inchis=convertCharArray(metsToAdd.inchis);
    if numel(metsToAdd.inchis)~=nMets
        EM='metsToAdd.inchis must have the same number of elements as metsToAdd.mets';
        dispEM(EM);
    end
    %Add empty field if it doesn't exist
    if ~isfield(newModel,'inchis')
        newModel.inchis=largeFiller;
    end
    newModel.inchis=[newModel.inchis;metsToAdd.inchis(:)];
else
    %Add empty strings if structure is in model
    if isfield(newModel,'inchis')
        newModel.inchis=[newModel.inchis;filler];
    end
end

if isfield(metsToAdd,'metSmiles')
    metsToAdd.metSmiles=convertCharArray(metsToAdd.metSmiles);
    if numel(metsToAdd.metSmiles)~=nMets
        EM='metsToAdd.metSmiles must have the same number of elements as metsToAdd.mets';
        dispEM(EM);
    end
    %Add empty field if it doesn't exist
    if ~isfield(newModel,'metSmiles')
        newModel.metSmiles=largeFiller;
    end
    newModel.metSmiles=[newModel.metSmiles;metsToAdd.metSmiles(:)];
else
    %Add empty strings if structure is in model
    if isfield(newModel,'metSmiles')
        newModel.metSmiles=[newModel.metSmiles;filler];
    end
end

if isfield(metsToAdd,'metFormulas')
    metsToAdd.metFormulas=convertCharArray(metsToAdd.metFormulas);
    if numel(metsToAdd.metFormulas)~=nMets
        EM='metsToAdd.metFormulas must have the same number of elements as metsToAdd.mets';
        dispEM(EM);
    end
    %Add empty field if it doesn't exist
    if ~isfield(newModel,'metFormulas')
        newModel.metFormulas=largeFiller;
    end
    newModel.metFormulas=[newModel.metFormulas;metsToAdd.metFormulas(:)];
else
    %Add default
    if isfield(newModel,'metFormulas')
        newModel.metFormulas=[newModel.metFormulas;filler];
    end
end

if isfield(metsToAdd,'metCharges')
    if numel(metsToAdd.metCharges)~=nMets
        EM='metsToAdd.metCharges must have the same number of elements as metsToAdd.mets';
        dispEM(EM);
    end
    if ~isnumeric(metsToAdd.metCharges)
        EM='metsToAdd.metCharges must be of type "double"';
        dispEM(EM);
    end
    if ~isfield(newModel,'metCharges')
        newModel.metCharges=NaN(numel(largeFiller),1);
    end
    newModel.metCharges=[newModel.metCharges;metsToAdd.metCharges(:)];
else
    %Add default
    if isfield(newModel,'metCharges')
        newModel.metCharges=[newModel.metCharges;NaN(numel(filler),1)];
    end
end

if isfield(metsToAdd,'metNotes')
    metsToAdd.metNotes=convertCharArray(metsToAdd.metNotes);
    if numel(metsToAdd.metNotes)==1 && numel(metsToAdd.mets)>1
        temp=cell(numel(metsToAdd.mets),1);
        temp(:)=metsToAdd.metNotes;
        metsToAdd.metNotes=temp;
    end
    if numel(metsToAdd.metNotes)~=nMets
        EM='metsToAdd.metNotes must have the same number of elements as metsToAdd.mets';
        dispEM(EM);
    end
    %Add empty field if it doesn't exist
    if ~isfield(newModel,'metNotes')
        newModel.metNotes=largeFiller;
    end
    newModel.metNotes=[newModel.metNotes;metsToAdd.metNotes(:)];
else
    %Add empty strings if structure is in model
    if isfield(newModel,'metNotes')
        newModel.metNotes=[newModel.metNotes;filler];
    end
end

%Don't check the type of metMiriams
if isfield(metsToAdd,'metMiriams')
    if numel(metsToAdd.metMiriams)==1 && numel(metsToAdd.mets)>1
        temp=cell(numel(metsToAdd.mets),1);
        temp(:)={metsToAdd.metMiriams};
        metsToAdd.metMiriams=temp;
    end
    if numel(metsToAdd.metMiriams)~=nMets
        EM='metsToAdd.metMiriams must have the same number of elements as metsToAdd.mets';
        dispEM(EM);
    end
    %Add empty field if it doesn't exist
    if ~isfield(newModel,'metMiriams')
        newModel.metMiriams=cell(nOldMets,1);
    end
    newModel.metMiriams=[newModel.metMiriams;metsToAdd.metMiriams(:)];
else
    if isfield(newModel,'metMiriams')
        newModel.metMiriams=[newModel.metMiriams;cell(nMets,1)];
    end
end

if isfield(newModel,'metFrom')
    newModel.metFrom=[newModel.metFrom;filler];
end

%Expand the S matrix
newModel.S=[newModel.S;sparse(nMets,size(newModel.S,2))];

if copyInfo==true
    [I, J]=ismember(metsToAdd.metNames,model.metNames);
    J=J(I);
    %I is the indexes of the new metabolites for which a metabolite with
    %the same name existed
    I=find(I)+nOldMets;
    %Go through each of the added mets and copy annotation if it doesn't
    %exist
    for i=1:numel(I)
        if isfield(newModel,'inchis')
            if isempty(newModel.inchis{I(i)})
                newModel.inchis(I(i))=newModel.inchis(J(i));
            end
        end
        if isfield(newModel,'metSmiles')
            if isempty(newModel.metSmiles{I(i)})
                newModel.metSmiles(I(i))=newModel.metSmiles(J(i));
            end
        end
        if isfield(newModel,'metFormulas')
            if isempty(newModel.metFormulas{I(i)})
                newModel.metFormulas(I(i))=newModel.metFormulas(J(i));
            end
        end
        if isfield(newModel,'metMiriams')
            if isempty(newModel.metMiriams{I(i)})
                newModel.metMiriams(I(i))=newModel.metMiriams(J(i));
            end
        end
        if isfield(newModel,'metCharges')
            newModel.metCharges(I(i))=newModel.metCharges(J(i));
        end
    end
end
end

author: SysBioChalmers

repo: RAVEN