text:
function msm_to_mm_test06 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST06 tests MSM_TO_MM_ARRAY_INTEGER_SKEW_SYMMETRIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST06\n' );
  fprintf ( 1, '  Convert an MSM to MM array integer skew-symmetric format.\n' );

  output_filename = 'msm_to_mm_test06.mm';

  a = i4mat_indicator ( 4, 4 );
  a = a - a';
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'array', 'integer', 'skew-symmetric' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function axisenlarge(f,h)
%AXISENLARGE - enlarge the axes of a figure (f,h)
%
% Usage:  (1) axisenlarge(1.05)    % enlarge axes by 5% in each direction
%         (2) axisenlarge(-1.05)   % shrink to fit content before enlarging
%
% Inputs:
%    f      enlarge axis by a factor f relative to current size or
%           by -f relative to the graph content. For separate factors
%           in each direction use [fx fy fz] or [fleft fright fbottom ftop fdown fup]
%    h      axis handle [default = gca]

%	   Copyright (C) Mike Brookes 2012
%      Version: $Id: axisenlarge.m 10428 2018-03-08 07:33:39Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fpt=[1 1 1 1 1 1; 1 1 2 2 2 2; 1 1 2 2 3 3; 1 2 3 4 3 4; 1 2 3 4 5 5; 1 2 3 4 5 6];
if nargin<2 || ~numel(h)
    h=gca;
end
if nargin<1 || ~numel(f)
    f=-1.02;
end
nf=min(numel(f),6);
f=f(fpt(nf,:));  % expand f to dimension 4
if any(f>=0)
    ax0=[get(h,'XLim') get(h,'YLim') get(h,'ZLim')];
else
    ax0=zeros(1,6);
end
if any(f<0)
    axis(h,'tight');
    ax1=[get(h,'XLim') get(h,'YLim')  get(h,'ZLim')];
    ax0(f<0)=ax1(f<0);
    f=abs(f);
end
ax1=ax0.*f+ax0([2 1 4 3 6 5]).*(1-f);
set(h,'XLim',ax1(1:2),'YLim',ax1(3:4),'ZLim',ax1(5:6));

author: jtkim-kaist

repo: Speech-enhancement


text:
%
% Multiply polynomials, polynomial matrices and numeric matrices
%             h = f1 * f2 * ... * fk
%
% Syntax:   (ptimes is the shortened alias for PolynomialTimes)
%            >> p = PolynomialTimes(f,g)
%            >> p = PolynomialTimes(f, g, h)  %
%   Input:    f, g, h, ...  --- (string/numeric/cell/matrix)
%                 polynomials, numeric matrices and/or polynomial cells,
%                 so that the multiplication is well-defined.
%
%  Output:    p --- (string) product of the input polynomials
%
%  Example:
%
%     >> PolynomialTimes(3,[1 5; -2 4],'x+y',{'y-1';'x+z'})
%
%     ans =
%
%    '-3*x + 15*x^2 - 3*y + 18*x*y + 3*y^2 + 15*x*z + 15*y*z'
%    '6*x + 12*x^2 + 6*y + 6*x*y - 6*y^2 + 12*x*z + 12*y*z'
%

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function CMap = cmap_cluster(nbVals)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015

CMap =[
    0.5569    0.7922    1.0000
    0.5419    0.7747    0.9804
    0.5270    0.7572    0.9608
    0.5121    0.7398    0.9412
    0.4972    0.7223    0.9216
    0.4822    0.7049    0.9020
    0.4673    0.6874    0.8824
    0.4524    0.6700    0.8627
    0.4374    0.6525    0.8431
    0.4225    0.6350    0.8235
    0.4076    0.6176    0.8039
    0.3927    0.6001    0.7843
    0.3777    0.5827    0.7647
    0.3628    0.5652    0.7451
    0.3479    0.5478    0.7255
    0.3330    0.5303    0.7059
    0.3180    0.5128    0.6863
    0.3031    0.4954    0.6667
    0.2882    0.4779    0.6471
    0.2732    0.4605    0.6275
    0.2583    0.4430    0.6078
    0.2434    0.4256    0.5882
    0.2285    0.4081    0.5686
    0.2135    0.3906    0.5490
    0.1986    0.3732    0.5294
    0.1837    0.3557    0.5098
    0.1688    0.3383    0.4902
    0.1538    0.3208    0.4706
    0.1389    0.3034    0.4510
    0.1240    0.2859    0.4314
    0.1090    0.2684    0.4118
    0.0941    0.2510    0.3922
    0.3922         0         0
    0.4118    0.0037    0.0023
    0.4314    0.0073    0.0046
    0.4510    0.0110    0.0068
    0.4706    0.0147    0.0091
    0.4902    0.0183    0.0114
    0.5098    0.0220    0.0137
    0.5294    0.0257    0.0159
    0.5490    0.0293    0.0182
    0.5686    0.0330    0.0205
    0.5882    0.0367    0.0228
    0.6078    0.0404    0.0250
    0.6275    0.0440    0.0273
    0.6471    0.0477    0.0296
    0.6667    0.0514    0.0319
    0.6863    0.0550    0.0342
    0.7059    0.0587    0.0364
    0.7255    0.0624    0.0387
    0.7451    0.0660    0.0410
    0.7647    0.0697    0.0433
    0.7843    0.0734    0.0455
    0.8039    0.0770    0.0478
    0.8235    0.0807    0.0501
    0.8431    0.0844    0.0524
    0.8627    0.0880    0.0546
    0.8824    0.0917    0.0569
    0.9020    0.0954    0.0592
    0.9216    0.0991    0.0615
    0.9412    0.1027    0.0638
    0.9608    0.1064    0.0660
    0.9804    0.1101    0.0683
    1.0000    0.1137    0.0706];
end

author: brainstorm-tools


text:
%%train network config
param.caffe_path=fullfile(fileparts(pwd),'matcaffe');
param.solver_netfile=fullfile(fileparts(pwd),'model','solver_bn.prototxt');
param.fintune_model=fullfile(fileparts(pwd),'pretrain_model','googlenet_bn.caffemodel');
param.test_net_file=fullfile(fileparts(pwd),'model','test_bn.prototxt');
param.test_batch_size=16;
param.result_save_file='result.txt';
param.save_model_file='model_cross';
param.save_model_name='base_train_cross';
param.test_interval=500;
param.train_maxiter=4000;
param.data_path=fullfile(fileparts(pwd),'generate_data');
param.traindata_filename='train_data_Prid';
param.testdata_filename='test_data_Prid';
param.use_data_split_index=1;
param.data_split_num=1;
param.train_person_num=150;
param.use_gpu=1;
param.gpu_id=0;
%%
for split_index=param.use_data_split_index:param.data_split_num
    if ~exist(strcat(param.save_model_file,num2str(split_index)),'file')
        mkdir(strcat(param.save_model_file,num2str(split_index)));
    end
    load(strcat(param.data_path,param.traindata_filename,num2str(split_index),'/train_data.mat'));
    %% find caffe
    cur_path=pwd;
    addpath(genpath(param.caffe_path))
    cd(param.caffe_path);
    caffe.reset_all;
    caffe.init_log(fullfile(cur_path,'log'));
    addpath(param.caffe_path);
    if param.use_gpu
        caffe.set_mode_gpu;
    else
        caffe.set_mode_cpu;
    end
    cd(cur_path);
    caffe_solver=caffe.get_solver(param.solver_netfile,param.gpu_id);
    net=caffe.get_net(param.test_net_file,'test');
    caffe_solver.use_caffemodel(param.fintune_model);
    input_data_shape=caffe_solver.nets{1}.blobs('data').shape;
    batch_size=input_data_shape(4);
    batch_data=zeros(input_data_shape,'single');
    batch_label=zeros(1,1,1,batch_size,'single');
    train_x_axis=[];
    train_y_axis=[];
    iter=0;
    while iter<param.train_maxiter

        [batch_data,batch_label]=get_train_minibatch(batch_data,batch_label,batch_size,train_data_cam1,train_data_cam2,label_train_cam1,label_train_cam2);
        caffe_solver.nets{1}.blobs('data').set_data(batch_data);
        caffe_solver.nets{1}.blobs('label').set_data(batch_label);
        caffe_solver.step(1);
        iter = caffe_solver.iter;
        if mod(iter,10)==0
            soft_loss=caffe_solver.nets{1}.blobs('loss3/loss').get_data;
            acc=caffe_solver.nets{1}.blobs('loss3/top-1').get_data;
            loss1=caffe_solver.nets{1}.blobs('loss1/loss').get_data;
            loss2=caffe_solver.nets{1}.blobs('loss2/loss').get_data;
%             loss4=caffe_solver.net.blobs('loss4').get_data;
            train_x_axis=[train_x_axis,iter];
            train_y_axis=[train_y_axis,soft_loss];
            plot(train_x_axis,train_y_axis);
            drawnow;
            fprintf('iter= %d ,softmaxloss= %f ,loss1=%f  loss2=%f\n',iter,soft_loss,loss1,loss2);
        end
        if mod(iter,param.test_interval)==0
            caffe_solver.nets{1}.save(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
            net.copy_from(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
            net.blobs('data').reshape([224 224 3 param.test_batch_size]);
            load(strcat(param.data_path,param.testdata_filename,num2str(split_index),'/test_data.mat'));
            test_script;
        end
        if iter>100 && mod(iter,100)==0
            caffe_solver.nets{1}.save(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
        end
    end
%     test_network;
end

author: liuyuisanai

repo: Quality-Aware-Network


text:
function zd=zdres_sat(flag,r,obs,nav,azel,dant,opt,rtk,zd) %#ok

global glc
lam=nav.lam(obs.sat,:); nf=rtk.NF;

if opt.ionoopt==glc.IONOOPT_IFLC
    if lam(1)==0||lam(2)==0,return;end

    % check snr mask,not surpport

    f1=glc.CLIGHT/lam(1); f2=glc.CLIGHT/lam(2);
    C1=f1^2/(f1^2-f2^2);  C2=-f2^2/(f1^2-f2^2);

    dant_if=C1*dant(1)+C2*dant(2);

    if obs.L(1)~=0&&obs.L(2)~=0
        zd.y(1)=C1*obs.L(1)*lam(1)+C2*obs.L(2)*lam(2)-r-dant_if;
    end
    if obs.P(1)~=0&&obs.P(2)~=0
        zd.y(2)=C1*obs.P(1)+C2*obs.P(2)-r-dant_if;
    end
else
    for i=1:nf
        if lam(i)==0,continue;end

        % check snr mask,not surpport

        if obs.L(i)~=0,zd.y(i   )=obs.L(i)*lam(i)-r-dant(i);end
        if obs.P(i)~=0,zd.y(i+nf)=obs.P(i)       -r-dant(i);end
    end
end

return

repo: GINav


text:
function out = read_nervus_data(nrvHdr, segment, range, chIdx)
% read_nervus_data  Returns data from Nicolet file.
%
%   OUT = read_nervus_data(NRVHDR, SEGMENT, RANGE, CHIDX) returns data in an n x m array of
%   doubles where n is the number of datapoints and m is the number
%   of channels.
%
%   NRVHDR is a header from the function read_nervus_header
%   SEGMENT is the segment number in the file to read from
%   RANGE is a 1x2 array with the [StartIndex EndIndex]  - default: all
%   and CHIDX is a vector of channel indeces - default: all
%
%   FILENAME is the file name of a file in the Natus/Nicolet/Nervus(TM)
%   format (originally designed by Taugagreining HF in Iceland)
%
%   Based on ieeg-portal/Nicolet-Reader
%   at https://github.com/ieeg-portal/Nicolet-Reader
%
% Copyright (C) 2016, Jan Brogger and Joost Wagenaar
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%  FieldTrip is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
%
%  FieldTrip is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: $

if nargin == 0
  ft_error('Missing argument');
elseif nargin == 1
  segment = 1;
  range = [1 nrvHdr.targetSampleCount];
  chIdx = nrvHdr.matchingChannels(:);
elseif nargin == 2
  range = [1 nrvHdr.targetSampleCount];
  chIdx = nrvHdr.matchingChannels(:);
elseif nargin == 3
  chIdx = nrvHdr.matchingChannels(:);
end

assert(length(range) == 2, 'Range is [firstIndex lastIndex]');
assert(range(1) > 0, 'Range is must start at one or higher');
assert(length(segment) == 1, 'Segment must be single value.');

% Get cumulative sum segments.
cSumSegments = [0 cumsum([nrvHdr.Segments.duration])];

% Reopen .e file.
h = fopen_or_error(nrvHdr.filename,'r','ieee-le');

% Find sectionID for channels
lChIdx = length(chIdx);
sectionIdx = zeros(lChIdx,1);
for i = 1:lChIdx
  tmp = find(strcmp(num2str(chIdx(i)-1),{nrvHdr.StaticPackets.tag}),1);
  sectionIdx(i) = nrvHdr.StaticPackets(tmp).index;
end

% Iterate over all requested channels and populate array.
out = zeros(range(2) - range(1) + 1, lChIdx);
for i = 1 : lChIdx
  %disp(['Reading channel ' num2str(chIdx) ' segment ' num2str(segment)]);
  % Get sampling rate for current channel
  curSF = nrvHdr.Segments(segment).samplingRate(chIdx(i));
  mult = nrvHdr.Segments(segment).scale(chIdx(i));

  % Find all sections
  allSectionIdx = nrvHdr.allIndexIDs == sectionIdx(i);
  allSections = find(allSectionIdx);

  % Find relevant sections
  sectionLengths = [nrvHdr.MainIndex(allSections).sectionL]./2;
  cSectionLengths = [0 cumsum(sectionLengths)];

  skipValues = cSumSegments(segment) * curSF;
  firstSectionForSegment = find(cSectionLengths > skipValues, 1) - 1 ;
  lastSectionForSegment = [];
  %lastSectionForSegment = firstSectionForSegment + ...
  %  find(cSectionLengths > curSF*nrvHdr.Segments(segment).duration,1) - 2 ;

  if isempty(lastSectionForSegment)
    lastSectionForSegment = length(cSectionLengths);
  end

  offsetSectionLengths = cSectionLengths - cSectionLengths(firstSectionForSegment);

  firstSection = find(offsetSectionLengths < range(1) ,1,'last');

  samplesInChannel = nrvHdr.Segments(segment).sampleCount;
  if range(2) > samplesInChannel
    endRange = samplesInChannel;
  else
    endRange = range(2);
  end

  lastSection = find(offsetSectionLengths >= endRange,1)-1;

  if isempty(lastSection)
    lastSection = length(offsetSectionLengths)-1;
  end

  if lastSection > lastSectionForSegment
    ft_error('Index out of range for current section: %i > %i, on channel: %i', ...
      range(2), cSectionLengths(lastSectionForSegment+1), chIdx(i));
  end

  useSections = allSections(firstSection: lastSection);
  useSectionL = sectionLengths(firstSection: lastSection);

  % First Partial Segment
  curIdx = 1;
  curSec = nrvHdr.MainIndex(useSections(1));
  fseek(h, curSec.offset,'bof');

  firstOffset = range(1) - offsetSectionLengths(firstSection);
  lastOffset = min([range(2) useSectionL(1)]);
  lsec = lastOffset-firstOffset + 1;

  fseek(h, (firstOffset-1) * 2,'cof');
  out(1 : lsec,i) = fread(h, lsec, 'int16') * mult;
  curIdx = curIdx +  lsec;

  if length(useSections) > 1
    % Full Segments
    for j = 2: (length(useSections)-1)
      curSec = nrvHdr.MainIndex(useSections(j));
      fseek(h, curSec.offset,'bof');

      out(curIdx : (curIdx + useSectionL(j) - 1),i) = ...
        fread(h, useSectionL(j), 'int16') * mult;
      curIdx = curIdx +  useSectionL(j);
    end

    % Final Partial Segment
    curSec = nrvHdr.MainIndex(useSections(end));
    fseek(h, curSec.offset,'bof');
    %Sometimes there is less data available in the file.
    %So we calculate the wanted lastReadLength,
    %but assign in the data array the actualReadLength
    lastReadLength = length(out)-curIdx+1;
    %disp(['LastReadLength: ' num2str(lastReadLength)]);
    lastData = fread(h, lastReadLength, 'int16');
    actualReadLength = length(lastData);
    %if (actualReadLength <lastReadLength)
    %  warning('Some data not found');
    %disp(['ActualLastReadLength: ' num2str(lastReadLength)]);
    out(curIdx : (curIdx+actualReadLength-1),i) = lastData * mult;
  end

end

% Close the .e file.
fclose(h);

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function b = det(a)
% function b=det(a)
%
% DESCRIPTION
%   Determinant of a matrix polynomial
%
% INPUTS
%   A: polynomial
%
% OUTPUTS
%   B: polynomial, the determinant of A
%
% SYNTAX
%   B = det(A);

% 6/14/2002: PJS  Initial Coding

sza=size(a);
if sza(1)~=sza(2)
    error('Matrix must be square');
end

if isempty(a)
    b = polynomial(1);
elseif sza(1)==1
    b = a;
else
    L.type = '()';
    b = polynomial(0);
    for i1 = 1:sza(1);
        % Recursive cofactor expansion for det
        % XXX Faster algos for computing det exist
        L.subs = {[2:sza(1)] [1:i1-1, i1+1:sza(1)]};
        M = subsref(a,L);
        cofactor = (-1)^(1+i1)*det(M);

        L.subs = {1, i1};
        a1_i1 = subsref(a,L);
        b = b + a1_i1*cofactor;
    end
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [candidates, score] = run_bing(im, num_candidates)
  bing_path = '/path/to/bing.sh';
  model_prefix = '/path/to/BING-Objectness/VOC2007/Results/ObjNessB2W8MAXBGR';

  [ret_code, tmp_dir] = system('mktemp -d');
  assert(ret_code == 0);
  tmp_dir = strtrim(tmp_dir);
  fprintf('created %s\n', tmp_dir);

  img_filename = fullfile(tmp_dir, 'foo.png');
  imwrite(im, img_filename);

  candidates_file = fullfile(tmp_dir, 'candidates.txt');

  cmd = sprintf('%s -i %s -o %s -m %s -d %d', ...
    bing_path, img_filename, candidates_file, model_prefix, num_candidates);
  [ret_code] = system(cmd);
  assert(ret_code == 0);

  [score, x1, y1, x2, y2] = textread(candidates_file, '%f %d %d %d %d', 'headerlines', 1);
  candidates = [x1 y1 x2 y2];

  system(sprintf('rm -rf "%s"', tmp_dir));
end

author: hosang

repo: detection-proposals


text:
% SOLVE_BILAPLACE_2D_ISO: Solve a 2d bilaplace problem with a variational formulation based on the Laplacian.
%
% The function solves the bilaplacian problem
%
%   laplace( epsilon(x) laplace(u)) = f    in Omega = F((0,1)^2)
%                             du/dn = 0    on Gamma_D
%                                 u = 0    on Gamma_D
%
% with a variational formulation based on the Laplacian, that is,
%
%       (epsilon laplace(u), laplace(v)) = (f,v)
%
% USAGE:
%
%  [geometry, msh, space, u] = solve_bilaplace_2d_iso (problem_data, method_data)
%
% INPUT:
%
%  problem_data: a structure with data of the problem. It contains the fields:
%    - geo_name:     name of the file containing the geometry
%    - drchlt_sides: sides with Dirichlet boundary condition (always homogeneous)
%    - c_diff:       physical parameter (epsilon in the equation)
%    - f:            source term
%
%  method_data : a structure with discretization data. Its fields are:
%    - degree:     degree of the spline functions.
%    - regularity: continuity of the spline functions.
%    - nsub:       number of subelements with respect to the geometry mesh
%                   (nsub=1 leaves the mesh unchanged)
%    - nquad:      number of points for Gaussian quadrature rule
%
% OUTPUT:
%
%  geometry: geometry structure (see geo_load)
%  msh:      mesh object that defines the quadrature rule (see msh_cartesian)
%  space:    space object that defines the discrete space (see sp_scalar)
%  u:        the computed degrees of freedom
%
% Copyright (C) 2009, 2010, 2011 Carlo de Falco
% Copyright (C) 2011, 2013, 2015 Rafael Vazquez
% Copyright (C) 2013, Marco Pingaro
%
%    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/>.
function [geometry, msh, space, u] = ...
              solve_bilaplace_2d_iso (problem_data, method_data)

% Extract the fields from the data structures into local variables
data_names = fieldnames (problem_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= problem_data.(data_names{iopt});']);
end
data_names = fieldnames (method_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= method_data.(data_names{iopt});']);
end

if (any (regularity < 1))
  error ('The regularity should be at least 1')
end

% Construct geometry structure
geometry = geo_load (geo_name);
degelev  = max (degree - (geometry.nurbs.order-1), 0);
nurbs    = nrbdegelev (geometry.nurbs, degelev);
[rknots, zeta, nknots] = kntrefine (nurbs.knots, nsub-1, nurbs.order-1, regularity);

nurbs = nrbkntins (nurbs, nknots);
geometry = geo_load (nurbs);

% Construct msh structure
rule     = msh_gauss_nodes (nquad);
[qn, qw] = msh_set_quad_nodes (zeta, rule);
msh      = msh_cartesian (zeta, qn, qw, geometry,'der2', true);

% Construct space structure
space = sp_nurbs (geometry.nurbs, msh);

% Assemble the matrices
stiff_mat = op_laplaceu_laplacev_tp (space, space, msh, c_diff);
rhs       = op_f_v_tp (space, msh, f);

% Apply boundary conditions
% Only homogeneous conditions are implemented
u = zeros (space.ndof, 1);

drchlt_dofs_u = []; drchlt_dofs_r = [];
for iside = drchlt_sides
  drchlt_dofs_u = union (drchlt_dofs_u, space.boundary(iside).dofs);
  drchlt_dofs_r = union (drchlt_dofs_r, space.boundary(iside).adjacent_dofs);
end
drchlt_dofs = union (drchlt_dofs_u, drchlt_dofs_r);

int_dofs = setdiff (1:space.ndof, drchlt_dofs);
rhs(int_dofs) = rhs(int_dofs) - stiff_mat(int_dofs, drchlt_dofs)*u(drchlt_dofs);

% Solve the linear system
u(int_dofs) = stiff_mat(int_dofs, int_dofs) \ rhs(int_dofs);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [ vX ] = SolveProxTvAdmm( vY, mD, paramLambda, numIterations )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveProxTvAdmm( vY, mD, paramLambda, numIterations )
% Solves the Prox of the Total Variation (TV) Norm using Alternating
% Direction Method of Multipliers (ADMM) Method.
% Basically solves the problem given by:
% $$ \arg \min_{ x \in \mathbb{R}^{n} } \frac{1}{2} {\left\| x - y \right|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$
% Input:
%   - vX                -   input Vector.
%                           Initialization of the iterative process.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - mA                -   Input Matirx.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L1 Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  Wikipedia ADMM - https://en.wikipedia.org/wiki/Augmented_Lagrangian_method#Alternating_direction_method_of_multipliers.
% Remarks:
%   1.  Using vanilla ADMM with no optimization of the parameter or
%       smoothing.
%   2.  Matrix Factorization caching according to "Matrix Inversion Lemma"
%       (See S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein -
%       Distributed Optimization and Statistical Learning via the
%       Alternating Direction Method of Multipliers Page 28). Basically:
%       (mA.' * mA + paramRho * I)^(-1) = (1 / paramRho) + (1 / (paramRho * paramRho)) * mA.' * (I + (1 /
%       paramRho) * mA * mA.')^(-1) * mA
% Known Issues:
%   1.  A
% TODO:
%   1.  Pre calculate decomposition of the Linear System.
% Release Notes:
%   -   1.0.000     27/11/2019  Royi Avital
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

paramRho = 5;

mI = eye(size(vY, 1));
mC = decomposition(mI + paramRho * (mD.' * mD), 'chol');

vX = vY;
vZ = ProxL1(mD * vX, paramLambda / paramRho);
vU = mD * vX - vZ;

for ii = 2:numIterations

    vX = mC \ (vY + (paramRho * mD.' * (vZ - vU)));
    vZ = ProxL1(mD * vX + vU, paramLambda / paramRho);
    vU = vU + mD * vX - vZ;

end

end

function [ vX ] = ProxL1( vX, lambdaFactor )

% Soft Thresholding
vX = max(vX - lambdaFactor, 0) + min(vX + lambdaFactor, 0);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function determ = cheby_van2_determinant ( n )

%*****************************************************************************80
%
%% CHEBY_VAN2_DETERMINANT returns the determinant of the CHEBY_VAN2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  if ( n <= 0 )
    determ = 0.0;
  elseif ( n == 1 )
    determ = 1.0;
  else
    determ = r8_mop ( floor ( n / 2 ) ) * sqrt ( 2.0 )^( 4 - n );
  end

  return
end


text:
function [soln,eqn,info] = PoissonWG(node,elem,bdFlag,pde,option)
%% POISSONWG Poisson equation: lowest order weak Galerkin element
%
%   u = POISSONWG(node,elem,bdFlag,pde) produces the linear finite element
%   approximation of the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N,
%       Robin boundary condition     g_R*u + d*grad(u)*n=g_N on \Gamma _R
%
% The usage is the same as Poisson. Weak Galerkin method on a triangle is
% summarized in PoissonWGfemrate for detail.
%
%   Example
%
%     squarePoissonWG;
%
%   See also Poisson, squarePoissonWG
%
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Preprocess
if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('option','var'), option = []; end
% important constants
NT = size(elem,1);
N = size(node,1);

%% Diffusion coefficient
time = cputime;  % record assembling time
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'dquadorder'), option.dquadorder = 1; end
if ~isempty(pde.d) && isnumeric(pde.d)
   K = pde.d;                                 % d is an array
end
if ~isempty(pde.d) && ~isnumeric(pde.d)       % d is a function
    [lambda,weight] = quadpts(option.dquadorder);
    nQuad = size(lambda,1);
    K = zeros(NT,1);
    for p = 1:nQuad
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
        K = K + weight(p)*pde.d(pxy);
   end
end

%% Construct data structure
[elem2edge,edge] = dofedge(elem);
NE = size(edge,1);
Ndof = NT + NE;
elem2dof = NT + elem2edge;

%% Assemble stiffness matrix
A = sparse(Ndof,Ndof);

% compute ct2 = 1/mean(||x-xc||^2)
center = (node(elem(:,1),:) + node(elem(:,2),:) + node(elem(:,3),:))/3;
mid1 = (node(elem(:,2),:) + node(elem(:,3),:))/2;
mid2 = (node(elem(:,3),:) + node(elem(:,1),:))/2;
mid3 = (node(elem(:,1),:) + node(elem(:,2),:))/2;
ct2 = 3./sum((mid1 - center).^2 + (mid2 - center).^2 + (mid3 - center).^2,2);
[Dphi,area] = gradbasis(node,elem);
clear center mid1 mid2 mid3

% Mbb: edge - edge
for i = 1:3
    for j = i:3
        % local to global index map
        ii = double(elem2dof(:,i));
        jj = double(elem2dof(:,j));
        % local stiffness matrix
        Aij = 4*dot(Dphi(:,:,i),Dphi(:,:,j),2).*area + 4/9*ct2.*area;
        if ~isempty(pde.d)
            Aij = K.*Aij;
        end
        if (j==i)
            A = A + sparse(ii,jj,Aij,Ndof,Ndof);
        else
            A = A + sparse([ii,jj],[jj,ii],[Aij; Aij],Ndof,Ndof);
        end
    end
end

% Mob: interior - edge
Aij = -4/3*ct2.*area;
if ~isempty(pde.d)
    Aij = K.*Aij;
end
Mob = sparse([(1:NT)', (1:NT)', (1:NT)'], ...
             double(elem2dof(:)), [Aij, Aij, Aij], Ndof, Ndof);
A = A + Mob + Mob';

% Moo: diagonal of interor
Aij = 4*ct2.*area;
if ~isempty(pde.d)
    Aij = K.*Aij;
end
A =  A + sparse(1:NT, 1:NT, Aij, Ndof, Ndof);

clear K Aij

%% Assemble the right hand side
b = zeros(Ndof,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 2;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if isreal(pde.f) % f is a real number or vector and not a function
   switch length(pde.f)
       case NT  % f is piecewise constant
         b(1:NT) = pde.f.*area;
       case N   % f is piecewise linear
         b(1:NT) = (pde.f(elem(:,1)) + pde.f(elem(:,2)) + pde.f(elem(:,3)))/3.*area;
       case 1   % f is a scalar e.g. f = 1
         b(1:NT) = pde.f*area;
   end
end
if ~isempty(pde.f) && ~isreal(pde.f)  % f is a function
    [lambda,weight] = quadpts(option.fquadorder);
	nQuad = size(lambda,1);
    bt = zeros(NT,1);
    for p = 1:nQuad
		% quadrature points in the x-y coordinate
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
		fp = pde.f(pxy);
        bt = bt + weight(p)*fp;
    end
    bt = bt.*area;
    b(1:NT) = bt;
end
clear pxy bt

%% Set up boundary conditions
if nargin<=3, bdFlag = []; end
[AD,b,u,freeDof,isPureNeumann] = getbdWG(A,b);

%% Record assembling time
assembleTime = cputime - time;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(freeDof), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if NE <= 1e3  % Direct solver for small size systems
        option.solver = 'direct';
    else          % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeDof) = AD(freeDof,freeDof)\b(freeDof);
        residual = norm(b - AD*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
        % eleminate elementwise dof
        option.solver = 'CG';
        if isfield(option,'reducesystem') && (option.reducesystem == 0)
            option.x0 = u;
            [u,info] = mg(AD,b,elem,option,edge);
        else
            option.x0 = u(NT+1:end);
            Aoinv = spdiags(1./diag(AD(1:NT,1:NT)),0,NT,NT);
            Aob = AD(1:NT,NT+1:end);
            Abo = Aob';
            Abb = AD(NT+1:end,NT+1:end);
            Abbm = Abb - Abo*Aoinv*Aob;
            bm = -Abo*Aoinv*b(1:NT) + b(NT+1:end);
            [ub,info] = mg(Abbm,bm,elem,option,edge);
            u(1:NT) = Aoinv*(b(1:NT) - Aob*ub);
            u(NT+1:end) = ub;
        end
    case 'amg'
        option.solver = 'CG';
        [u(freeDof),info] = amg(AD(freeDof,freeDof),b(freeDof),option);
end
% post-process for pure Neumann problem
if isPureNeumann
    uc = sum(u(1:NT).*area);
    u = u - uc;   % normalization for pure Neumann problem
end

%% Compute Du
dudx =  u(elem2dof(:,1)).*Dphi(:,1,1) + u(elem2dof(:,2)).*Dphi(:,1,2) ...
      + u(elem2dof(:,3)).*Dphi(:,1,3);
dudy =  u(elem2dof(:,1)).*Dphi(:,2,1) + u(elem2dof(:,2)).*Dphi(:,2,2) ...
      + u(elem2dof(:,3)).*Dphi(:,2,3);
Du = -2*[dudx, dudy];

%% Output information
if nargout == 1
    soln = u;
else
    soln = struct('u',u,'Du',Du);
    eqn = struct('A',AD,'b',b,'edge',edge,'freeDof',freeDof);
    info.assembleTime = assembleTime;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdWG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [AD,b,u,freeDof,isPureNeumann]= getbdWG(A,b)
    %% GETBDCR Boundary conditions for Poisson equation: WG element.

    u =zeros(Ndof,1);
    %% Initial check
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    idxR = (bdFlag(:) == 3);      % index of Robin edges in bdFlag
    if any(idxR)
        isRobin = false(NE,1);
        isRobin(elem2edge(idxR)) = true;
        Robin = edge(isRobin,:);  % Robin edges
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        ve = node(Robin(:,1),:) - node(Robin(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        mid = (node(Robin(:,1),:) + node(Robin(:,2),:))/2;
        ii = NT + find(isRobin);  % for WG: edge dof is after elem dof
        ss = pde.g_R(mid).*edgeLength; % exact for linear g_R
        A = A + sparse(ii,ii,ss,Ndof,Ndof);
    end

    % Find Dirichlet boundary nodes: fixedEdge
    fixedEdge = []; freeEdge = [];
    if ~isempty(bdFlag)              % find boundary edges
        idxD = (bdFlag(:) == 1);     % all Dirichlet edges in bdFlag
        isFixedEdge = false(NE,1);
        isFixedEdge(elem2edge(idxD)) = true;  % index of fixed boundary edges
        fixedEdge = find(isFixedEdge);
        freeEdge = find(~isFixedEdge);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        % no bdFlag, only pde.g_D is given in the input
        s = accumarray(elem2edge(:), 1, [NE 1]);
        fixedEdge = find(s == 1);
        freeEdge = find(s == 2);
    end
    isPureNeumann = false;
    if isempty(fixedEdge) && isempty(Robin)  % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        isPureNeumann = true;
        fixedEdge = 1;
        freeEdge = (2:NE)';    % eliminate the kernel by enforcing u(1) = 0;
    end
    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedEdge,fixedEdge)=I, AD(fixedEdge,freeEdge)=0, AD(freeEdge,fixedEdge)=0.
    if ~isempty(fixedEdge)
        bdidx = zeros(Ndof,1);
        bdidx(NT + fixedEdge) = 1;
        Tbd = spdiags(bdidx,0,Ndof,Ndof);
        T = spdiags(1-bdidx,0,Ndof,Ndof);
        AD = T*A*T + Tbd;
    else
        AD = A;
    end

    %% Part 2: Find boundary edges and modify the right hand side b
    % Find boundary edges: Neumann
    Neumann = [];
    if ~isempty(bdFlag)  % bdFlag specifies different bd conditions
        idxN = (bdFlag(:) == 2);      % all Neumann edges in bdFlag
        isNeumann = elem2edge(idxN | idxR); % index of Neumann and Robin edges
        % since boundary integral is also needed for Robin edges
        Neumann = edge(isNeumann,:);      % Neumann edges
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        s = accumarray(elem2edge(:), 1, [NE 1]);
        Neumann = edge(s == 1,:);
    end

    % Neumann boundary condition
    if ~isempty(Neumann) && ~isempty(pde.g_N) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 3;   % default order exact for linear gN
        end
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        ge = zeros(size(Neumann,1),1);
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        for pp = 1:nQuadgN
            % quadrature points in the x-y coordinate
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNp = pde.g_N(ppxy);
            ge = ge+ weightgN(pp)*gNp;
        end
        ge = ge.*el;
        b(NT+isNeumann) = b(NT+isNeumann) + ge;
    end
    % The case with non-empty Neumann edges but g_N=0 or g_N=[] corresponds to
    % the zero flux boundary condition on Neumann edges and no modification of
    % A,u,b is needed.

    % Dirichlet boundary condition
    if ~isPureNeumann && ~isempty(fixedEdge) && ...
       ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && all(pde.g_D == 0))    % nonzero g_D
        if isnumeric(pde.g_D)  % pde.g_D could be a numerical array
            u(fixedEdge) = pde.g_D(fixedEdge);
        else % pde.g_D is a function handle
            mid = (node(edge(fixedEdge,1),:) + node(edge(fixedEdge,2),:))/2;
            u(NT+fixedEdge) = pde.g_D(mid);
        end
        b = b - A*u;
        b(NT+fixedEdge) = u(NT+fixedEdge);
    end
    % The case with non-empty Dirichlet nodes but g_D=0 or g_D=[] corresponds
    % to the zero Dirichlet boundary condition and no modification of u,b is
    % needed.

    % Pure Neumann boundary condition
    if isPureNeumann
        b = b - mean(b);   % compatilbe condition: sum(b) = 0
        b(1) = 0;
    end

    freeDof = [(1:NT)'; NT+freeEdge];
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % end of PoissonWG

repo: ifem


text:
function handle = draw_gmm(haxes,Priors,Mu,Sigma )
%DRAW_GMM Summary of this function goes here
%   Detailed explanation goes here

K = size(Priors,2);
STD=3;

for k=1:K
    l = max(Priors);
    for i=3:(-1):STD
        w = Priors(k)/l;
        handle = plot_gaussian_ellipsoid(Mu(:,k),Sigma(:,:,k),i,100,haxes,w,[0 0 1]);
        set(handle,'LineWidth',2);
      %  obj.text_handle = [obj.text_handle text(obj.Mu(1,k),obj.Mu(2,k),num2str(k),'FontWeight','bold','FontSize',24,'Color',obj.color,'VerticalAlignment','middle','HorizontalAlignment','center')];
    end
end

end

author: epfl-lasa


text:
% This files adds a coastline from an existing data set
global coastline faults mainfault main
report_this_filefun(mfilename('fullpath'));

%aa = a;
[file1,path1] = uigetfile( '*.mat',' Earthquake Datafile'); %disabled window position
loadpath = [path1 file1];

new_data = load(loadpath);

loaded=false;
if isfield(new_data,'coastline')
    if ~isempty(new_data.coastline)
        loaded=true;
        coastline=new_data.coastline;
    end
end

if isfield(new_data,'faults')
    if ~isempty(new_data.faults)
        loaded=true;
        faults=new_data.faults;
    end
end

if isfield(new_data,'mainfault')
    if ~isempty(new_data.mainfault)
        loaded=true;
        mainfault=new_data.mainfault;
    end
end

if isfield(new_data,'main')
    if ~isempty(new_data.main)
        loaded=true;
        main=new_data.main;
    end
end

if ~loaded
    disp('Error lodaing data! Are they in the right *.mat format??')
end
whos
%a = aa;
%clear aa
update(mainmap())

author: CelsoReyes


text:
function [ results ] = run_DAT_AFFINE( seq, res_path, bSaveImage )

project_dir = fileparts(fileparts(fileparts(mfilename('fullpath'))));
addpath(fullfile(project_dir, 'trackers', 'DAT', 'src'));

visualization = 0;
use_normal_size = 1;
normal_height = 360;
normal_width = 640;
I = imread(seq.s_frames{1});
if use_normal_size
  scale_height = size(I, 1) / normal_height;
  scale_width = size(I, 2) / normal_width;
else
  scale_height = 1;
  scale_width = 1;
end
% Load annotations
groundtruth = seq.init_rect;
groundtruth = groundtruth ./ [scale_width, scale_height, scale_width, scale_height];

% Load default settings
cfg = affine_parameters_dat();

% Tracking
frames = 1:length(seq.s_frames);
times = zeros(size(frames));
rects = zeros(length(frames), 4);
spf = tic;
for frame = frames
  I = imread(seq.s_frames{frame});
  % I = imresize(I, [360, 640]);
  if use_normal_size
    I = imresize(I, [normal_height, normal_width]);
  end

  ttrack = tic;
  if frame == 1
    [state, location] = tracker_dat_initialize(I, groundtruth, cfg);
  else
    [state, location] = tracker_dat_update(state, I, cfg);
  end
  times(frame) = toc(ttrack);
  % Visualization
  % rects(frame, :) = location*2;
  rects(frame, :) = location .* [scale_width, scale_height, scale_width, scale_height];
  if visualization
    figure(1), clf
    imshow(I)
    rectangle('Position', location, 'EdgeColor', 'b', 'LineWidth', 2);
    drawnow
  end
end
spf = toc(spf);

results.res = rects;
results.type='rect';
results.fps = length(frames)/spf;

end

author: flyers

repo: drone-tracking


text:
mex tamEStepInnerMex.c 3rd-party/util.c;
cd utils;
mex termCountsToWordListMex.c;
cd ..;
cd 3rd-party/fastfit
install_fastfit;
cd ../lightspeed
install_lightspeed;
cd ../..

author: jacobeisenstein

repo: SAGE


text:
function x = from_array(A,opt)
    %FROM_ARRAY Approximate full array by TTeMPS tensor of prescribed rank
    %   or within a prescribed tolerance.
    %
    %   X = TTeMPS.from_array( A, tol ) approximates the given array A by a
    %       TTeMPS tensor such that the the error is in the order of tol.
    %
    %   X = TTeMPS.from_array( A, r ), with r a vector of length (ndims(A))+1),
    %       approximates the given array A by a rank-r TTeMPS tensor, such that
    %       X.rank = r.
    %

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    n = size(A);
    d = length(n);

    if length(opt) == 1
        useTol = true;
        tol = opt;
        r = ones(1,d+1);
    else
        useTol = false;
        r = opt;
        if r(1) ~= 1 || r(d+1) ~= 1
            error('Invalid rank specified')
        end
    end

    U = cell(1,d);

    % process from left to right
    % first core
    A = reshape( A, n(1), prod(n(2:end)));
    [u,s,v] = svd(A,'econ');
    if useTol
        r(2) = trunc_singular( diag(s), tol );
    end
    U{1} = reshape( u(:,1:r(2)), [1, n(1), r(2)] );
    A = s(1:r(2),1:r(2))*v(:,1:r(2))';

    % middle cores
    for i = 2:d-1
        A = reshape( A, n(i)*r(i), prod(n(i+1:end)));
        [u,s,v] = svd(A,'econ');
        if useTol
            r(i+1) = trunc_singular( diag(s), tol );
        end
        U{i} = reshape( u(:,1:r(i+1)), [r(i), n(i), r(i+1)] );
        A = s(1:r(i+1),1:r(i+1)) * v(:,1:r(i+1))';
    end

    %last core
    U{d} = reshape(A, [r(d), n(d), 1]);

    x = TTeMPS( U );
end

author: NicolasBoumal

repo: manopt


text:
%GetNextItem - Get next item from a Batch object (iterator mode).
%
% This is a helper class to easily read and parse batch files.
%
%  USAGE
%
%    [b,item,line] = GetNextItem(b);
%
%    b              batch object

% Copyright (C) 2007 by Michaël Zugaro
%
% 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.

function [b,item,line] = GetNextItem(b)

% Check number of parameters
if nargin ~= 1,
	error('Incorrect number of parameters (type ''help GetNextItem'' for details).');
end

b.currentItem = b.currentItem+1;
b.currentField = 0;
if b.currentItem > size(b.field,1),
	item = [];
	line = [];
else
	item = b.currentItem;
	line = b.line(item);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function returnedmappings = getDatabaseMappings(field,qualifiers)
% getDataBaseMappings returns information on known mappings of database
% entries to model field names, along with additional information about the
% fields.
%
% INPUT:
%
%    field:                 the basic model field to extract mappings for (e.g. 'met', 'gene', 'rxn')
%
% OPTIONAL INPUT:
%
%    qualifiers:            the qualifiers to restrict the selection to.
%                           These have to be part of the bioql modifiers
%                           definition (e.g. is, isDescribedBy, isEncodedBy
%                           etc) providing 'all' will return all associated
%                           Database mappings. Default: 'all'
%
% OUTPUT:
%
%    returnedmappings:      The mappings known for the given field. The
%                           structure is:
%                           X{:,1} : the database ID (in identifiers.org/miriam annotation)
%                           X{:,2} : The qualifier associated with the DB
%                           X{:,3} : The model field associated with this db
%                           X{:,4} : The association field (met/rxn/gene/prot/comp)
%                           X{:,5} : The specified regular expression the identifier has to adhere to.
%                           X{:,6} : The type of the qualifier (modelQualifier or bioQualifier)

if exist('qualifiers', 'var') && ...
    ((ischar(qualifiers) && isequal(qualifiers,'all')) || iscell(qualifiers) && any(ismember(qualifiers,'all')))
    %If it exists and either is 'all' or contains 'all', return all.
    qualifiers = 'all';
elseif ~exist('qualifiers','var')
    %if it doesn't exist also return all
    qualifiers = 'all';
end

mappings = getDefinedFieldProperties('Database',true);

returnedmappings = mappings(strcmp(field,mappings(:,4)),:);
%If not all, restrict
if ~isequal(qualifiers,'all')
    returnedmappings = returnedmappings(ismember(returnedmappings(:,2),qualifiers),:);
end

author: opencobra

repo: cobratoolbox


text:
function [mu, s2, MM, VV] = simulLMGPmc(hyp, inffcn, meanfcn, covfcn, likfcn, input, target, targetvariance,...
    derivinput, derivtarget, derivvariance, xt, lag, Nsamples)
% simulLMGPmc - Simulation of the dynamic GP model with incorporated local models (LMGP models),
% where the output variance is propagated using the numerical Monte Carlo approximation
%
%% Syntax
%  [mu, s2, MM, VV] = simulLMGPmc(hyp, inf, mean, cov, lik, input, target, targetvariance,...
%     derivinput, derivtarget, derivvariance, xt, lag, Nsamples)
%
%% Description
% Idea: at every time step the output of GP model is approximated with
% the Gaussian distribution, from which we sample one value. We repeat the
% procedure Nsamples-times.
% Nsamples samples, which are used as the future inputs of the GP model.
% Samples are re-used if necessary (ie. y(k-1) for y(k-2) if lag=2 etc.)
% Currently it can be used only with the Gaussian covariance function and
% with the white noise model (sum of covSEard and covNoise) due to the gpSD00.
% Uses routine gpSD00.
% see K. Azman. Identifikacija dinamicnih sistemov z Gaussovimi procesi. PhD
% thesis, University of Ljubljana, Ljubljana, 2007 (in Slovene).
%
% Input:
% * hyp            ... the structure of hyperparameters
% * inf      	   ... the inference method 	      --> not used, for interface compatibility only
% * cov      	   ... the prior covariance function  --> not used, for interface compatibility only
% * mean    	   ... the prior mean function        --> not used, for interface compatibility only
% * lik      	   ... the likelihood function        --> not used, for interface compatibility only
% * input          ... the input part of the training data,  NxD matrix
% * target         ... the output part of the training data (ie. target), Nx1 vector
% * targetvariance ... the target variance, use NaN where not known
% * derivinput     ... the input part of the derivative training data, NEQxD matrix
% * derivtarget    ... target derivatives, NEQxD matrix
% * derivvariance  ... variances of the local model prameters, NEQxD matrix
% * xt             ... the input matrix for simulation, kxD vector, see
%                      construct_ARXsimul_input.m for more info
% * lag            ... the order of the model (number of used lagged outputs)
% * Nsamples       ... the number of samples used in algorithm (ie. runs of simulation)
%
% Output:
% * mu             ... the mean predicted output
% * s2             ... associated variances (including noise variance)
% * MM             ... the matrix of all predicted means, kxNsamples
% * VV             ... associated predicted variances
%
% See Also
% gpSD00.m, simulLMGPnaive.m, simulGPmc.m
%
% Examples
% demo_example_lmgp_simulation.m
%
%%
% * Written by K.Azman, 31.05.2005
% * Based on the work of C.E. Rasmussen and A. Girard.
%
%
% Changelog:
%
% 16.2.2015, Martin Stepancic:
%		 	-changed the function interface as gpml > 3.0
%			-removed the addition of autocovariance - this is now
%			 already included in gpSD00.m.
%

fun_name = 'simulLMGPmc';

[n, D]      = size(input);
fullinput  = [input; repmat(derivinput,D,1)];
fulltarget = [target; derivtarget(:)];
[N, D]      = size(fullinput);
[nD, D] = size(derivinput);

MM = zeros(Nsamples,size(xt,1));
VV = zeros(Nsamples,size(xt,1));

for jjj=1:Nsamples

    if(mod(jjj,10)==0)
        disp([fun_name, ', run:  ',int2str(jjj),'/',int2str(Nsamples)]);
    end

    % 1st point - input is "point"
    test = xt(1,:);
    [mu(1), s2(1)] = gpSD00(hyp, inffcn, meanfcn, covfcn, likfcn, input, target, ...
			targetvariance, derivinput, derivtarget, derivvariance, test);

    for k=2:length(xt)

        test = [xt(k, lag+1:end)];

        % For the NEXT prediction...
        % assumed normal distribution, for more accurate procedure see
        % simulGPmc

        % random sample from assumed distribution
        ysampled = mu(k-1) + randn(1)*sqrt(s2(k-1));
        % generate input for the GP model
        if (k>lag)
            test = [mu(k-lag:k-2) ysampled xt(k, lag+1:end)];
        elseif (k<=lag)
            test = [xt(k, 1:lag-k+1) mu(1:k-2) ysampled xt(k, lag+1:end)];
        end

        [mu(k), s2(k)] = gpSD00(hyp, inffcn, meanfcn, covfcn, likfcn, input, target, ...
			targetvariance, derivinput, derivtarget, derivvariance, test);
    end

    MM(jjj,:) = mu;
    VV(jjj,:) = s2;

end
% individual realisations saved in matrices MM and VV
% approximate all output distributions with Gaussian distribution
mu = mean(MM);
s2 = mean(VV) + mean((MM-repmat(mu,Nsamples,1)).^2);

mu = mu';
s2 = s2';

MM = MM';
VV = VV';

return

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function [Bj1p,PSp] = FBASE( B,a)
config;
fS = nnconfig.FilterSize;
s=fS*fS-1;
size = nnconfig.ImageSize ;

if nargin == 1
    for j = 1:s
        Bj = reshape(B(:,j),fS,fS);
        Bj1(:,:,j) = rot90(Bj,2);
        PS(:,:,j) = psf2otf(Bj1(:,:,j), size);
    end
    PSp =gpuArray(PS);
    Bj1p = gpuArray(Bj1);
end
if nargin == 2
    for j = 1:s
        Bj = reshape(B(:,j),fS,fS);
        Bj1(:,:,j) = rot90(Bj,2);
    end
    Bj1p = gpuArray(Bj1)  ;
end
end

repo: Deep-ADMM-Net


text:
function HDR=openeep(HDR,arg2,arg3,arg4,arg5,arg6)
% OPENEEP opens EEProbe files (but does not read the data).
% However, it is recommended to use SOPEN instead of OPENEEP.
%
% see also: SLOAD, SOPEN, SREAD, SCLOSE, SEOF, STELL, SSEEK.

%	$Id: openeep.m 2205 2009-10-27 12:18:15Z schloegl $
%	Copyright (c) 2007 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the  License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

if strcmp(HDR.TYPE,'EEProbe-CNT'),

       	HDR.FILE.FID = fopen(HDR.FileName,'rb');
	H = openiff(HDR.FILE.FID);

	if isfield(H,'RIFF');
                HDR.RIFF      = H.RIFF;
                HDR.FILE.OPEN = 1;
                HDR.Label     = {};
                HDR.PhysDim   = {};
                HDR.SPR       = inf;
                HDR.NRec      = 1;
                HDR.FILE.POS  = 0;
                if ~isfield(HDR.RIFF,'CNT');
                	HDR.TYPE = 'unknown';
                elseif ~isfield(HDR.RIFF.CNT,'eeph') | ~isfield(HDR.RIFF.CNT,'LIST');
                	HDR.TYPE = 'unknown';
                else
			s = char(HDR.RIFF.CNT.eeph);
			field = '';
			while ~isempty(s)
				[line,s] = strtok(s,[10,13]);
				if strncmp(line,'[Sampling Rate]',15);
					field = 'SampleRate';
				elseif strncmp(line,'[Samples]',9);
					field = 'SPR';
				elseif strncmp(line,'[Channels]',10);
					field = 'NS';
				elseif strncmp(line,'[Basic Channel Data]',20);
					k = 0;
					while (k<HDR.NS),
						[line,s] = strtok(s,[10,13]);
						if ~strncmp(line,';',1);
							k = k+1;
							[num,status,sa]=str2double(line);
							HDR.Label{k} = sa{1};
							HDR.PhysDim{k,1} = sa{4};
							HDR.Cal(k) = num(2)*num(3);
						end;
					end;
				elseif strncmp(line,'[History]',9);
					[t,s] = strtok(s,[10,13]);
					while ~strncmp(t,'EOH',3)
						[t,s] = strtok(s,[10,13]);
					end;
				elseif strncmp(line,';',1);
				elseif strncmp(line,'[',1);
					field = '';
				elseif ~isempty(field);
					[num,status,sa] = str2double(line);
					if ~status,
						HDR = setfield(HDR,field,num);
						field = '';
					end;
				end;
			end;

			% decode data block
			HDR.SPR = H.RIFF.CNT.LIST.raw3.ep(1);
			HDR.NRec = length(H.RIFF.CNT.LIST.raw3.ep)-2;

if 0, %try
	%%% ### FIXME ### %%%
	%%% decoding of data block
	%%% this is work in progress.
			HDR.EEP.epoch_length = H.RIFF.CNT.LIST.raw3.ep(1);
			HDR.EEP.epoch_start  = [H.RIFF.CNT.LIST.raw3.ep(2:end),length(H.RIFF.CNT.LIST.raw3.data)];
			n = [HDR.EEP.epoch_length,length(H.RIFF.CNT.LIST.raw3.chan)];
%dec2bin(double(HDR.RIFF.CNT.LIST.raw3.data(HDR.EEP.epoch_start(1:end-1)+1)),8)
			%n=n([2,1]);
			accu  = zeros(1,HDR.NS);
			for k = 1:HDR.NRec,
				ix1 = HDR.EEP.epoch_start(k)+1;
				bytes = double(H.RIFF.CNT.LIST.raw3.data(HDR.EEP.epoch_start(k)+1:HDR.EEP.epoch_start(k+1)));
				ix1 = 1;
			for ch= 1:HDR.NS,
				ix1 = ceil(ix1);
				byte1 = H.RIFF.CNT.LIST.raw3.data(HDR.EEP.epoch_start(k)+ix1);
				meth  = double(bitshift(byte1,-4));

				if any(meth==[1:3])
					nbits = bitand(bytes(1),15);
					nexcbits = bitshift(bytes(2),4);
					if ~nexcbits, nexcbits = 16; end;
					y = repmat(NaN,HDR.SPR,1);
					y(1) = bitshift(bitand(bytes(2),15)*256+bytes(3),4)+bitshift(bytes(4),-4);
					bitpos = 28;
					ix1 = ix1 + 3.5;
					for k1 = 2:HDR.SPR,
						%ix1 = bitpos/8;
						ix2 = ix1 + nbits/8;
						i1 = ceil(ix1);
						r1 = (ix1-i1+1)*8;
						i2 = ceil(ix2);
						r2 = (ix2-i2+1)*8;
%[r1,i1,r2,i2],
						a = bitand(bytes(i1),2^r1-1);
						a = bitshift(a,16) + bitshift(bytes(i1+1),8) + bytes(i1+2);
						a = bitshift(a,r2-16-r1);
						ix1 = ix2;
						if a>2^(nbits-1)
							a=a-2^nbits;
						elseif a==2^(nbits-1)
							ix2 = ix1 + nexcbits/8;
							i1 = ceil(ix1);
							r1 = (ix1-i1+1)*8;
							i2 = ceil(ix2);
							r2 = (ix2-i2+1)*8;
							a = bitand(bytes(i1),2^r1-1);
							a = bitshift(a,16) + bitshift(bytes(i1+1),8) + bytes(i1+2);
							a = bitshift(a,r2-16-r1);
							ix1 = ix2;
							if a>2^(nexcbits-1)
								a=a-2^nexcbits;
							end;
						end;
						y(k1)=a;
					end;
					%fprintf(HDR.FILE.stderr,'Warning SOPEN(EEProbe): decompression method %i not implented\n',meth);
				elseif any(meth==[9:11])
					nbits = bitshift(bitand(bytes(1),15),2) + bitshift(bytes(2),-6);
					nexcbits = bitand(bytes(2),63);
					y = repmat(NaN,HDR.SPR,1);
					y(1) = bytes(3:6)*(256.^[3,2,1,0]');
					bitpos = 48;
					ix1 = ix1 + 6;
					for k1 = 2:HDR.SPR,
						%ix1 = bitpos/8;
						ix2 = ix1 + nbits/8;
						i1 = ceil(ix1);
						r1 = (ix1-i1+1)*8;
						i2 = ceil(ix2);
						r2 = (ix2-i2+1)*8;
%[r1,i1,r2,i2],
						a = bitand(bytes(i1),2^r1-1);
						a = bitshift(a,32) + bitshift(bytes(i1+1),24) + bitshift(bytes(i1+2),16) + bitshift(bytes(i1+3),8) + bytes(i1+4);
						a = bitshift(a,r2-32-r1);
						ix1 = ix2;
						if a>2^(nbits-1)
							a=a-2^nbits;
						elseif a==2^(nbits-1)
							ix2 = ix1 + nexcbits/8;
							i1 = ceil(ix1);
							r1 = (ix1-i1+1)*8;
							i2 = ceil(ix2);
							r2 = (ix2-i2+1)*8;
							a = bitand(bytes(i1),2^r1-1);
							a = bitshift(a,16) + bitshift(bytes(i1+1),8) + bytes(i1+2);
							a = bitshift(a,r2-16-r1);
							ix1 = ix2;
							if a>2^(nexcbits-1)
								a=a-2^nexcbits;
							end;
						end;
						y(k1)=a;
					end;
					fprintf(HDR.FILE.stderr,'Warning SOPEN(EEProbe): decompression method %i not implemented\n',meth);
					error('####');
				elseif meth==0,
					y = bytes(ix1+(1:2:HDR.SPR*2));
					y = y - 256*(y>127);
					y = 256*y + bytes(ix1+(2:2:HDR.SPR*2));
					ix1 = ix1 + 1 + 2 * HDR.SPR;
				elseif meth==8,
					y = bytes(ix1+(1:HDR.SPR*4));
					y(1:4:end) = y(1:4:end) - 256*(y(1:4:end)>127);
					y = reshape(y,[4,length(y)/4])'*(256.^[3;2;1;0]);
					ix1 = ix1 + 1 + 4 * HDR.SPR;
				else
					y = repmat(NaN,HDR.SPR,1);
					%fprintf(HDR.FILE.stderr,'ERROR SOPEN(EEProbe): decompression method %i not supported\n',double(meth));
					error(sprintf('EEProbe: decompression method %d not supported\n',double(meth)));
				end;

				if any(meth==[1,9]);
					y = reshape(y,[HDR.SPR,1]);
					y(1,:) = accu(ch)+y(1,:);
					y = cumsum(y,1);
					accu(ch) = y(end,:);
				elseif any(meth==[2,10]);
					y = reshape(y,[HDR.SPR,1]);
					y(2,:) = accu(ch)+y(1,:);
					y(2:end,:) = cumsum(y(2:end,:),1);
					accu(ch) = y(end,:);
				elseif any(meth==[3,11]);
					y = cumsum(reshape(y,[HDR.SPR,1]),2);
				elseif any(meth==[0,8]);
					%y = reshape(y,n([2,1]))';
					y = reshape(y,[HDR.SPR,1]);
				end;
				HDR.data(HDR.SPR*(k-1)+1:HDR.SPR*k,H.RIFF.CNT.LIST.raw3.chan(ch)+1) = y;
			end;
			end;
			HDR.Calib = sparse(2:HDR.NS+1,1:HDR.NS,HDR.Cal);
			HDR.TYPE = 'native';
else %catch
			HDR.Calib = sparse(2:HDR.NS+1,1:HDR.NS,1); % because SREAD uses READ_EEP_CNT.MEX
end;
                end
                HDR.SPR = HDR.SPR*HDR.NRec;
                HDR.NRec = 1;
        end
end;

	% read event file, if applicable
	fid = 0;
	if strcmp(HDR.TYPE,'EEProbe-TRG'),
	        fid = fopen(HDR.FileName,'rt');
        elseif strcmp(HDR.TYPE,'EEProbe-CNT')
	        fid = fopen(fullfile(HDR.FILE.Path,[HDR.FILE.Name,'.trg']),'rt');
	end;
	if fid>0,
                tmp = str2double(fgetl(fid));
                if ~isnan(tmp(1))
                	HDR.EVENT.SampleRate = 1/tmp(1);
	                N = 0;
	                while ~feof(fid),
	                        tmp = fscanf(fid, '%f %d %s', 3);
	                        if ~isempty(tmp)
	                                N = N + 1;
	                                HDR.EVENT.POS(N,1)  = round(tmp(1)*HDR.EVENT.SampleRate);
	                                HDR.EVENT.TYP(N,1)  = 0;
	                                %HDR.EVENT.DUR(N,1) = 0;
	                                %HDR.EVENT.CHN(N,1) = 0;

	                                HDR.EVENT.TeegType{N,1} = char(tmp(3:end));
	                                HDR.EVENT.TYP(N,1)  = str2double(HDR.EVENT.TeegType{N,1});		% numeric
	                        end
	                end;
	                HDR.EVENT.TYP(isnan(HDR.EVENT.TYP))=0;
	                HDR.TRIG = HDR.EVENT.POS(HDR.EVENT.TYP>0);
	%                HDR.EVENT.POS = HDR.EVENT.POS(HDR.EVENT.TYP>0);
	%                HDR.EVENT.TYP = HDR.EVENT.TYP(HDR.EVENT.TYP>0);
	%                HDR.EVENT = rmfield(HDR.EVENT,'TeegType');
		end;
                fclose(fid);
        end;

	if strcmp(HDR.TYPE,'EEProbe-AVR'),
	        % it appears to be a EEProbe file with an averaged ERP
        	try
        	        tmp = read_eep_avr(HDR.FileName);
        	catch
        	        fprintf(HDR.FILE.stderr,'ERROR SOPEN (EEProbe): Cannot open EEProbe-file, because read_eep_avr.mex not installed. \n');
        	        fprintf(HDR.FILE.stderr,'ERROR SOPEN (EEProbe): see http://www.smi.auc.dk/~roberto/eeprobe/\n');
        	        return;
	        end

        	% convert the header information to BIOSIG standards
	        HDR.FILE.FID = 1;               % ?
        	HDR.FILE.POS = 0;
	        HDR.NS = tmp.nchan;             % number of channels
        	HDR.SampleRate = tmp.rate;      % sampling rate
	        HDR.NRec  = 1;                   % it is an averaged ERP, therefore one record
        	HDR.SPR   = tmp.npnt;             % total number of samples in the file
	        HDR.Dur   = tmp.npnt/tmp.rate;    % total duration in seconds
        	HDR.Calib = [zeros(1,HDR.NS) ; eye(HDR.NS, HDR.NS)];  % is this correct?
	        HDR.Label = char(tmp.label);
	        HDR.PhysDim   = 'uV';
	        HDR.FLAG.UCAL = 1;
	        HDR.FILE.POS  = 0;
	        HDR.AS.endpos = HDR.SPR;
	        HDR.Label = tmp.label;
	        HDR.TriggerOffset = 0;

	        HDR.EEP.data = tmp.data';
	end;

author: PatternRecognition

repo: OpenBMI


text:
function [filtered_signal,filtb,filta]=lopass_butterworth(inputsignal,cutoff_freq,Fs,order)
% Low-pass Butterworth filter
% [filtered_signal,filtb,filta] = lopass_butterworth(inputsignal,cutoff_freq,Fs,order)
%
% This is simply a set of built-in Matlab functions, repackaged for ease of
% use by Chad Greene, October 2012.
%
% INPUTS:
% inputsignal = input time series
% cutoff_freq = filter corner frequency
% Fs = data sampling frequency
% order = order of Butterworth filter
%
% OUTPUTS:
% filtered_signal = the filtered time series
% filtb, filta = filter numerator and denominator (optional)
%
% EXAMPLE 1:
% load train
% t = (1:length(y))/Fs;
% y_filt = lopass_butterworth(y,900,Fs,4); % cut off at 900 Hz
% figure
% plot(t,y,'b',t,y_filt,'r')
% xlabel('time in seconds')
% box off
% legend('unfiltered','filtered')
% sound(y,Fs)      % play original time series
% pause(2)         % pause two seconds
% sound(y_filt,Fs) % play filtered time series
%
%
% EXAMPLE 2:
% load train
% t = (1:length(y))/Fs;
% [y_filt,filtb,filta] = lopass_butterworth(y,900,Fs,4); % cut off at 900 Hz
% [h1,f1] = freqz(filtb,filta,256,Fs);
%
% figure
% subplot(3,1,1)
% plot(t,y,'b',t,y_filt,'r')
% xlabel('time in seconds')
% box off
% text(0,.1,' time series','units','normalized')
%
% subplot(3,1,2)
% AX = plotyy(f1,10*log10(abs(h1)),f1,angle(h1),'semilogx');
% set(get(AX(1),'ylabel'),'string','gain (dB)')
% set(get(AX(2),'ylabel'),'string','phase (rad)')
% xlim(AX(1),[min(f1) max(f1)])
% xlim(AX(2),[min(f1) max(f1)])
% text(0,.1,' filter response','units','normalized')
% box off
%
% [Pxx,f] = pwelch(y,512,256,[],Fs,'onesided');
% [Pxxf,f_f]= pwelch(y_filt,512,256,[],Fs,'onesided');
% subplot(3,1,3)
% semilogx(f,10*log10(Pxx))
% hold on
% semilogx(f_f,10*log10(Pxxf),'r')
% xlabel('frequency (Hz)')
% ylabel('PSD (dB)')
% xlim([min(f1) max(f1)])
% box off
% legend('unfiltered','filtered','location','northwest')
% legend boxoff

nyquist_freq = Fs/2;  % Nyquist frequency
Wn=cutoff_freq/nyquist_freq;    % non-dimensional frequency
[filtb,filta]=butter(order,Wn,'low'); % construct the filter
filtered_signal=filtfilt(filtb,filta,inputsignal); % filter the data with zero phase

author: Sable

repo: mcbench-benchmarks


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.
%
% Note: This was slightly modified such that it expects y = 1 or y = 0

% Find Indices of Positive and Negative Examples
pos = find(y == 1); neg = find(y == 0);

% Plot Examples
plot(X(pos, 1), X(pos, 2), 'k+','LineWidth', 1, 'MarkerSize', 7)
hold on;
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7)
hold off;

end

author: Ayatans

repo: Machine-Learning-homework


text:
function y=voicebox(f,v)
%VOICEBOX  set global parameters for Voicebox functions Y=(FIELD,VAL)
%
%  Inputs:  F   is a field name
%           V   is a new value for the field
%
% Outputs:  Y   is set equal to the structure of parameters if the
%               f and v inputs are both present or both absent. If only
%               input f is specified, then y is set to the value of the
%               corresponding field or null if it doesn't exist.
%
% This routine contains default values for constants that are used by
% other functions in the VOICEBOX toolbox. Values in the first section below,
% entitled "System-dependent directory paths" should be set as follows:
%    PP.dir_temp     directory for storing temporary files
%    PP.dir_data     default directory to preappend to speech data file names
%                    when the "d" option is specified in READWAV etc.
%    PP.shorten      location of SHORTEN executable. SHORTEN is a proprietary file compression
%                    algorithm that is used for some SPHERE-format files. READSPH
%                    will try to call an external decoder if it is asked to
%                    read such a compressed file.
%    PP.sfsbin       location of Speech Filing Sysytem binaries. If the "c" option
%                    is given to READSFS, it will try to create a requested item
%                    if it is not present in the SFS file. This parameter tells it
%                    where to find the SFS executables.
%    PP.sfssuffix    suffix for Speech Filing Sysytem binaries. READSFS uses this paremeter
%                    to create the name of an SFS executable (see PP.sfsbin above).
% Other values defined in this routine are the defaults for specific algorithm constants.
% If you want to change these, please refer to the individual routines for a fuller description.

% Bugs/Suggestions
%    (1)  Could allow a * at the end of F to act as a wildcard and return/print a part structure

%      Copyright (C) Mike Brookes 2003
%      Version: $Id: voicebox.m,v 1.17 2011/07/28 06:48:45 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

persistent PP
if isempty(PP)

    % System-dependent directory paths and constants

    PP.dir_temp='F:\TEMP';                      % directory for storing temporary files
    PP.dir_data='E:\dmb\data\speech';           % default directory to preappend to speech data file names
    PP.shorten='C:\bin\shorten.exe';            % location of shorten executable
    PP.flac='C:\bin\flac.exe';                  % location of flac executable
    PP.sfsbin='F:\Program Files\SFS\Program';   % location of Speech Filing Sysytem binaries
    PP.sfssuffix='.exe';                        % suffix for Speech Filing Sysytem binaries
    PP.memsize=50e6;                            % Maximum amount of temporary memory to use (Bytes)

    % DYPSA glottal closure identifier

    PP.dy_cpfrac=0.3;           % presumed closed phase fraction of larynx cycle
    PP.dy_cproj=0.2;            % cost of projected candidate
    PP.dy_cspurt=-0.45;         % cost of a talkspurt
    PP.dy_dopsp=1;              % Use phase slope projection (1) or not (0)?
    PP.dy_ewdly=0.0008;         % window delay for energy cost function term [~ energy peak delay from closure] (sec)
    PP.dy_ewlen=0.003;          % window length for energy cost function term (sec)
    PP.dy_ewtaper=0.001;        % taper length for energy cost function window (sec)
    PP.dy_fwlen=0.00045;        % window length used to smooth group delay (sec)
    PP.dy_fxmax=500;            % max larynx frequency (Hz)
    PP.dy_fxmin=50;             % min larynx frequency (Hz)
    PP.dy_fxminf=60;            % min larynx frequency (Hz) [used for Frobenius norm only]
    PP.dy_gwlen=0.0030;         % group delay evaluation window length (sec)
    PP.dy_lpcdur=0.020;         % lpc analysis frame length (sec)
    PP.dy_lpcn=2;               % lpc additional poles
    PP.dy_lpcnf=0.001;          % lpc poles per Hz (1/Hz)
    PP.dy_lpcstep=0.010;        % lpc analysis step (sec)
    PP.dy_nbest=5;              % Number of NBest paths to keep
    PP.dy_preemph=50;           % pre-emphasis filter frequency (Hz) (to avoid preemphasis, make this very large)
    PP.dy_spitch=0.2;           % scale factor for pitch deviation cost
    PP.dy_wener=0.3;            % DP energy weighting
    PP.dy_wpitch=0.5;           % DP pitch weighting
    PP.dy_wslope=0.1;           % DP group delay slope weighting
    PP.dy_wxcorr=0.8;           % DP cross correlation weighting
    PP.dy_xwlen=0.01;           % cross-correlation length for waveform similarity (sec)

    % RAPT pitch tracker

    PP.rapt_f0min=50;           % Min F0 (Hz)
    PP.rapt_f0max=500;          % Max F0 (Hz)
    PP.rapt_tframe=0.01;        % frame size (s)
    PP.rapt_tlpw=0.005;         % low pass filter window size (s)
    PP.rapt_tcorw=0.0075;       % correlation window size (s)
    PP.rapt_candtr=0.3;         % minimum peak in NCCF
    PP.rapt_lagwt=0.3;          % linear lag taper factor
    PP.rapt_freqwt=0.02;        % cost factor for F0 change
    PP.rapt_vtranc=0.005;       % fixed voice-state transition cost
    PP.rapt_vtrac=0.5;          % delta amplitude modulated transition cost
    PP.rapt_vtrsc=0.5;          % delta spectrum modulated transition cost
    PP.rapt_vobias=0.0;         % bias to encourage voiced hypotheses
    PP.rapt_doublec=0.35;       % cost of exact doubling or halving
    PP.rapt_absnoise=0;         % absolute rms noise level
    PP.rapt_relnoise=2;         % rms noise level relative to noise floor
    PP.rapt_signoise=0.001;     % ratio of peak signal rms to noise floor (0.001 = 60dB)
    PP.rapt_ncands=20;          % max hypotheses at each frame
    PP.rapt_trms=0.03;                      % window length for rms measurement
    PP.rapt_dtrms=0.02;                     % window spacing for rms measurement
    PP.rapt_preemph=-7000;                  % s-plane position of preemphasis zero
    PP.rapt_nfullag=7;                      % number of full lags to try (must be odd)

    % now check some of the key values for validity

    if exist(PP.dir_temp)~=7        % check that temp directory exists
        PP.dir_temp = winenvar('temp');     % else use windows temp directory
    end

    [fnp,fnn,fne]=fileparts(mfilename('fullpath'));
    if exist(PP.shorten)~=2        % check that shorten executable exists
        PP.shorten=fullfile(fnp,'shorten.exe'); % next try local directory
        if exist(PP.shorten)~=2        % check if it exists in local directory
            PP.shorten='shorten.exe'; % finally assume it is on the search path
        end
    end

    if exist(PP.flac)~=2        % check that flac executable exists
        PP.flac=fullfile(fnp,'flac.exe'); % next try local directory
        if exist(PP.flac)~=2        % check if it exists in local directory
            PP.shorten='flac.exe'; % finally assume it is on the search path
        end
    end

end
if nargin==0
    if nargout==0
        % list all fields
        nn=sort(fieldnames(PP));
        cnn=char(nn);
        fprintf('%d Voicebox parameters:\n',length(nn));

        for i=1:length(nn);
            if ischar(PP.(nn{i}))
                fmt='  %s = %s\n';
            else
                fmt='  %s = %g\n';
            end
            fprintf(fmt,cnn(i,:),PP.(nn{i}));
        end
    else
        y=PP;
    end
elseif nargin==1
    if isfield(PP,f)
        y=PP.(f);
    else
        y=[];
    end
else
    if isfield(PP,f)
        PP.(f)=v;
        y=PP;
    else
        error(sprintf('''%s'' is not a valid voicebox field name',f));
    end
end

author: decouples


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% finishes a test run, nice output, small break, close all
%==============================================================================

function testEnd

caller = dbstack; % identify the name of the calling function
caller = caller(min(length(caller),2)).name;
caller = which(caller);

str = {
  sprintf('%% %s\n',char(ones(1,78)*'#'))
  sprintf('    file <%s> terminated\n',caller)
  sprintf('%% %s\n',char(ones(1,78)*'#'))
};

fprintf('%s\n',str{:});

builtin('pause',2);
close all

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


text:
function gn = g01_size ( )

%*****************************************************************************80
%
%% G01_SIZE returns the size for grid 1.
%
%  Discussion:
%
%    The "grid size" is simply the number of points in the grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer GN, the grid size.
%
  gn = 100;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ x, more ] = sgmga_vcn ( n, w, x, q_min, q_max, more )

%*****************************************************************************80
%
%% SGMGA_VCN returns the next constrained vector.
%
%  Discussion:
%
%    This function is intended to replace the "naive" version, now called
%    SGMGA_VCN_NAIVE, which is too slow for high dimensional problems.
%
%    For nonnegative vectors X of dimension N, and nonnegative
%    weights W, we define:
%
%      Q = sum ( 1 <= I <= N ) W(I) * X(I)
%
%    and seek X satisfying the constraint:
%
%      Q_MIN < Q <= Q_MAX
%
%    This routine returns, one at a time exactly those X which satisfy
%    the constraint.  No attempt is made to return the X values in
%    any particular order as far as Q goes.
%
%  Example:
%
%        W               4.0 3.0 5.0
%      MIN     16.0       0   0   0
%      ---     ----      -----------
%        1     20.0       5   0   0
%        2     19.0       4   1   0
%        3     18.0       3   2   0
%        4     17.0       2   3   0
%        5     20.0       2   4   0
%        6     19.0       1   5   0
%        7     18.0       0   6   0
%        8     17.0       3   0   1
%        9     20.0       3   1   1
%       10     19.0       2   2   1
%       11     18.0       1   3   1
%       12     17.0       0   4   1
%       13     20.0       0   5   1
%       14     18.0       2   0   2
%       15     17.0       1   1   2
%       16     20.0       1   2   2
%       17     19.0       0   3   2
%       18     19.0       1   0   3
%       19     18.0       0   1   3
%       20     20.0       0   0   4
%      ---     ----      ----------
%      MAX     20.0       6   7   5
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    An Anisotropic Sparse Grid Stochastic Collocation Method for Partial
%    Differential Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2411-2442.
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, real W(N), the weights, which should be nonnegative.
%    At least one weight must be positive.
%
%    Input, integer X(DIM_NUM).  On first call, with
%    MORE = FALSE, the input value of X is not important.  On subsequent calls,
%    the input value of X should be the output value from the previous call.
%
%    Input, real Q_MIN, Q_MAX, the lower and upper
%    limits on the sum.
%
%    Input, logical MORE.  On input, if the user has set MORE
%    FALSE, the user is requesting the initiation of a new sequence
%    of values.  If MORE is TRUE, then the user is requesting "more"
%    values in the current sequence.
%
%    Output, integer X(DIM_NUM).
%    On output, (with MORE = TRUE), the value of X will be the "next"
%    vector in the reverse lexicographical list of vectors that satisfy
%    the condition.  However, on output with MORE = FALSE, the vector
%    X is meaningless, because there are no more vectors in the list.
%
%    Output, logical MORE.  If MORE is TRUE on output,
%    then another value was found and returned in X, but if MORE is
%    FALSE, then there are no more values in the sequence, and X is
%    NOT the next value.
%
  persistent dir;
  persistent n2;
  persistent nstart;
  persistent xmax;
  persistent xmin;
%
%  Initialization for first call.
%
%  Allocate XMAX to remember the currently maximum possible value for each X.
%
%  Locate NSTART, the index of the first nonzero weight.
%  The algorithm is easier to program if the last index we look at
%  has a nonzero weight, so that it can always make up the remainder.
%
  if ( ~more )

    xmax = zeros ( n, 1 );
    nstart = - 1;

    for i = 1 : n
      if ( 0.0 < w(i) )
        nstart = i;
        break
      end
    end
%
%  Theoretically, we could even handle the case where all weights are zero.
%  That case is ruled out elsewhere in this software, so I will not try
%  to deal with it here for now.
%
    if ( nstart == - 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SGMGA_VCN - Fatal error!\n' );
      fprintf ( 1, '  No weight is positive.\n' );
      error ( 'SGMGA_VCN - Fatal error!' );
    end
%
%  Initialize X to zero, even the indices we ignore.
%
    x = zeros ( n, 1 );
%
%  N2 points to our current index of interest.
%
    n2 = n + 1;
    dir = - 1;

    more = 1;

  end
%
%  Look for the next solution vector X.
%
  while ( 1 )
%
%  If no more, the search is terminated.
%
    if ( ~more )

      break
%
%  DIR = -1, decrement N2, and, if possible, set X(N2) to XMIN.
%  DIR =  0, hold N2 at current value, and see if we can increment X(N2).
%
    elseif ( dir == - 1 || dir == 0 )

      if ( dir == - 1 )
        n2 = n2 - 1;
      end

      if ( w(n2) == 0.0 )

        xmin = 0;
        xmax(n2) = 0;

      elseif ( nstart < n2 )

        xmin = 0;
        xmax(n2) = floor ( ( q_max -  w(n2+1:n,1)' * x(n2+1:n,1) ) / w(n2) );

      elseif ( n2 == nstart && dir == - 1 )

        xmin = ceil ( ( q_min - w(n2+1:n,1)' * x(n2+1:n,1) ) / w(n2) );
        xmin = max ( xmin, 0 );
        if ( w(1:n2-1,1)' * x(1:n2-1,1) ...
           + w(n2) * xmin ...
           + w(n2+1:n,1)' * x(n2+1:n,1)  ...
          <= q_min )
          xmin = xmin + 1;
        end

        x(n2) = xmin;

        xmax(n2) = floor ( ( q_max -  w(n2+1:n,1)' * x(n2+1:n,1) ) / w(n2) );

      end

      if ( xmax(n2) < xmin )

        dir = + 1;

      else

        if ( n2 == nstart )

          if ( dir == - 1 )
            dir = 0;
            break
          elseif ( dir == 0 )
            x(n2) = x(n2) + 1;
            if ( x(n2) <= xmax(n2) )
              break
            else
              dir = + 1;
            end

          end

        else

          x(n2) = xmin;

        end

      end
%
%  DIR = + 1:
%  Try moving backwards to find an index N2 whose X we can increment.
%
    elseif ( dir == + 1 )

      while ( 1 )

        if ( n2 == n )
          dir = 0;
          more = 0;
          xmax = [];
          break;
        end

        n2 = n2 + 1;

        if ( 0.0 < w(n2) )

          if ( x(n2) < xmax(n2) )
            x(n2) = x(n2) + 1;
            dir = - 1;
            break
          end

        end
      end
    end
  end

  return
end


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

function [optval, optstd] = CallPut2(S,K,C)
% S = NSim x Nt matrix of simulated prices
% K = Strike price
% C = 1 -> Call; C = 0 -> Put

    Nt = size(S,2);         % Number of Discretization Steps
    V = S(:,Nt);
    if(C==1)
        optval = mean(max(V-K,0));
        optstd = std(max(V-K,0));
    else
        optval = mean(max(K-V,0));
        optstd = std(max(K-V,0));
    end
end


text:
% h = drawfly(trk,varargin)
% h = drawfly(x,y,theta,quartermaj,quartermin,varargin)
% draw triangle corresponding to fly position
% extra arguments will be fed to patch
function varargout = drawfly(x,y,varargin)

% draw an isosceles triangle with center (x,y)
% with rotation theta
% with height maj*4
% with base min*4

if isstruct(x),
  fly = x;
  t = y;
  x = fly.x(t);
  y = fly.y(t);
  theta = fly.theta(t);
  maj = fly.a(t);
  min = fly.b(t);
else
  if nargin < 5,
    error('not enough arguments; usage: drawflyo(x,y,theta,a,b,...)');
  end
  theta = varargin{1};
  maj = varargin{2};
  min = varargin{3};
  varargin = varargin(4:end);
end

if 0,

h = ellipsedraw(maj*2,min*2,x,y,theta);

else

% isosceles triangle not yet rotated or centered
pts = [-maj*2,-min*2,
       -maj*2,min*2,
       maj*2,0];

% rotate
costheta = cos(theta);
sintheta = sin(theta);
R = [costheta,sintheta;-sintheta,costheta];
pts = pts*R;

% translate
pts(:,1) = pts(:,1) + x;
pts(:,2) = pts(:,2) + y;

% plot
if isempty(varargin),
  h = patch(pts([1:3,1],1),pts([1:3,1],2),'b');
else
  h = patch(pts([1:3,1],1),pts([1:3,1],2),varargin{:});
end

end

if nargout > 0,
  varargout{1} = h;
end;

author: kristinbranson

repo: JAABA


text:
function X = cell2mat_ov_sum(I,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap,sz,Bs)

% converts a cell array to a matrix when the cell elements overlap
% INPUTS:
% I:            cell array
% grid_size:    true size of each element
% overlap:      amount of overlap in each direction
% d1:           number of rows of matrix
% d2:           number of columns of matrix

% OUTPUT:
% X:            output matrix

% Written by Eftychios A. Pnevmatikakis, Simons Foundation, 2016

if nargin < 10 || isempty(Bs)
    Bs = cellfun(@(x) ones(size(x)), I,'un',0);
end
X = zeros([sz,size(I{1,1},length(sz)+1)]);
B = zeros(size(X));
if length(sz) == 2; sz(3) = 1; end

for i = 1:length(xx_f)
    for j = 1:length(yy_f)
        for k = 1:length(zz_f)
            extended_grid = [max(xx_s(i)-overlap(1),1),min(xx_f(i)+overlap(1),sz(1)),max(yy_s(j)-overlap(2),1),min(yy_f(j)+overlap(2),sz(2)),max(zz_s(k)-overlap(3),1),min(zz_f(k)+overlap(3),sz(3))];
            %W = construct_weights([xx_s(i),xx_f(i),yy_s(j),yy_f(j),zz_s(k),zz_f(k)],extended_grid)';
            Xtemp = zeros(size(I{i,j,k}));
            Btemp = Xtemp;
            ind = ~isnan(I{i,j,k});
            Xtemp(ind) = Bs{i,j,k}(ind).*I{i,j,k}(ind);
            Btemp(ind) = Bs{i,j,k}(ind);
            %X(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6)) = X(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6)) + Bs{i,j,k}.*I{i,j,k};
            %B(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6)) = B(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6)) + Bs{i,j,k}.*(I{i,j,k}~=0);
            X(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6)) = X(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6)) + Xtemp;
            B(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6)) = B(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6)) + Btemp;
        end
    end
end

X = X./B;
%X(isnan(X))=0;

author: flatironinstitute

repo: NoRMCorre


text:
function x = false(x)
% FALSE Constrains a binary variable to be false (0)
%
% FALSE(x) returns the constraint x<=0.5.
%
%   See also SDPVAR/TRUE, SDPVAR/AND, SDPVAR/OR, SDPVAR/NOT, BINVAR, BINARY

x = (x<=0.5);

author: zarathustr

repo: LibQPEP


text:
function f = bohach2 ( m, x )

%*****************************************************************************80
%
%% BOHACH2 evaluates the Bohachevsky function #2.
%
%  Discussion:
%
%    The minimizer is
%
%      X* = [ 0.0, 0.0 ]
%      F(X*) = 0.0
%
%    Suggested starting point:
%
%      X = [ 0.6, 1.3 ];
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, integer M, the number of variables.
%
%    Input, real X(M), the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  f =       x(1) * x(1) ...
    + 2.0 * x(2) * x(2) ...
    - 0.3 * cos ( 3.0 * pi * x(1) ) * cos ( 4.0 * pi * x(2) ) ...
    + 0.3;

  return
end


text:
classdef VBox < uix.Box
    %uix.VBox  Vertical box
    %
    %  b = uix.VBox(p1,v1,p2,v2,...) constructs a vertical box and sets
    %  parameter p1 to value v1, etc.
    %
    %  A vertical box lays out contents from top to bottom.
    %
    %  See also: uix.HBox, uix.Grid, uix.VButtonBox, uix.VBoxFlex

    %  Copyright 2009-2016 The MathWorks, Inc.
    %  $Revision: 1594 $ $Date: 2018-03-28 02:27:52 +1100 (Wed, 28 Mar 2018) $

    properties( Access = public, Dependent, AbortSet )
        Heights % heights of contents, in pixels and/or weights
        MinimumHeights % minimum heights of contents, in pixels
    end

    properties( Access = protected )
        Heights_ = zeros( [0 1] ) % backing for Heights
        MinimumHeights_ = zeros( [0 1] ) % backing for MinimumHeights
    end

    methods

        function obj = VBox( varargin )
            %uix.VBox  Vertical box constructor
            %
            %  b = uix.VBox() constructs a horizontal box.
            %
            %  b = uix.VBox(p1,v1,p2,v2,...) sets parameter p1 to value v1,
            %  etc.

            % Set properties
            try
                uix.set( obj, varargin{:} )
            catch e
                delete( obj )
                e.throwAsCaller()
            end

        end % constructor

    end % structors

    methods

        function value = get.Heights( obj )

            value = obj.Heights_;

        end % get.Heights

        function set.Heights( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''Heights'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                ~any( isnan( value ) ), 'uix:InvalidPropertyValue', ...
                'Elements of property ''Heights'' must be real and finite.' )
            assert( isequal( size( value ), size( obj.Contents_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''Heights'' must match size of contents.' )

            % Set
            obj.Heights_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.Heights

        function value = get.MinimumHeights( obj )

            value = obj.MinimumHeights_;

        end % get.MinimumHeights

        function set.MinimumHeights( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''MinimumHeights'' must be of type double.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                all( value >= 0 ), 'uix:InvalidPropertyValue', ...
                'Elements of property ''MinimumHeights'' must be non-negative.' )
            assert( isequal( size( value ), size( obj.Heights_ ) ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''MinimumHeights'' must match size of contents.' )

            % Set
            obj.MinimumHeights_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.MinimumHeights

    end % accessors

    methods( Access = protected )

        function redraw( obj )
            %redraw  Redraw
            %
            %  c.redraw() redraws the container c.

            % Compute positions
            bounds = hgconvertunits( ancestor( obj, 'figure' ), ...
                [0 0 1 1], 'normalized', 'pixels', obj );
            heights = obj.Heights_;
            minimumHeights = obj.MinimumHeights_;
            padding = obj.Padding_;
            spacing = obj.Spacing_;
            r = numel( heights );
            xPositions = [padding + 1, max( bounds(3) - 2 * padding, 1 )];
            xPositions = repmat( xPositions, [r 1] );
            ySizes = uix.calcPixelSizes( bounds(4), heights, ...
                minimumHeights, padding, spacing );
            yPositions = [bounds(4) - cumsum( ySizes ) - padding - ...
                spacing * transpose( 0:r-1 ) + 1, ySizes];
            positions = [xPositions(:,1), yPositions(:,1), ...
                xPositions(:,2), yPositions(:,2)];

            % Set positions
            children = obj.Contents_;
            for ii = 1:numel( children )
                uix.setPosition( children(ii), positions(ii,:), 'pixels' )
            end

        end % redraw

        function addChild( obj, child )
            %addChild  Add child
            %
            %  c.addChild(d) adds the child d to the container c.

            % Add to sizes
            obj.Heights_(end+1,:) = -1;
            obj.MinimumHeights_(end+1,:) = 1;

            % Call superclass method
            addChild@uix.Box( obj, child )

        end % addChild

        function removeChild( obj, child )
            %removeChild  Remove child
            %
            %  c.removeChild(d) removes the child d from the container c.

            % Remove from sizes
            tf = obj.Contents_ == child;
            obj.Heights_(tf,:) = [];
            obj.MinimumHeights_(tf,:) = [];

            % Call superclass method
            removeChild@uix.Box( obj, child )

        end % removeChild

        function reorder( obj, indices )
            %reorder  Reorder contents
            %
            %  c.reorder(i) reorders the container contents using indices
            %  i, c.Contents = c.Contents(i).

            % Reorder
            obj.Heights_ = obj.Heights_(indices,:);
            obj.MinimumHeights_ = obj.MinimumHeights_(indices,:);

            % Call superclass method
            reorder@uix.Box( obj, indices )

        end % reorder

    end % template methods

end % classdef

author: goGPS-Project


text:
%% EVENT OBJECT (event.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This class is designed to define a generic agent and import this variables
% into the simulation space for the purpose of multi-vehicle control simulation.

% Author: James A. Douthwaite

classdef eventDefinition
%%% EVENT BASE CLASS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This class contains the basic properties of a generic event, neither
    properties
        % EVENT PARAMETERS
        eventID;
        name;
        type;
        time;
        info;
    end
%%  CLASS METHODS
    methods
%       CONSTRUCTION METHOD
       function obj = eventDefinition(time,name,summaryInfo)
            % This function is designed to generate an event with an
            % automatically defined ID number and basic properties
            % INPUTS:
            % time
            % name
            % summaryInfo - Assigned description string
            % OUTPUTS:
            % obj     - The generated object

            % INPUT HANDLING
            % CONFIRM EVENT TIME
            if ~isnumeric(time)                    % Default time setting
                %             if ~exist('time','var') || ~isnumeric(time)                    % Default time setting
                obj.time = 0;
                %warning('An event time must be defined.');
                return
            else
                obj.time = time;
            end
            % CONFIRM NAME STRING
            if isnumeric(name) || isempty(name)    % Default name setting
                obj.name = 'EVENT';
            else
                obj.name = name;
            end
            % CONFIRM DEFAULT SUMMARY INFORMATION
            %             if ~exist('summaryInfo','var') || isempty(summaryInfo)
            if isempty(summaryInfo)
                obj.info = 'None.';                                        % Default information setting
            else
                obj.info = summaryInfo;
            end

            % ALLOCATE AUTOINCREMENTING ID NUMBER
            persistent eventcount;
            if isempty(eventcount)
                eventcount = 1;
            else
                eventcount = eventcount + 1;
            end
            obj.eventID = uint16(eventcount);   % Assign the number as an ID tag

            % ASSIGN THE DEFAULT TYP
            obj.type = eventType.event;         % Assign default enumeration of 'event'
       end
    end
    % PRIVATE METHODS (CLASS SPECIFIC TOOLS)
    methods (Access = private)
        % PRINT THE OBJECT SUMMARY
        function displaySummary(obj)
            fprintf('[%s]\tevent created (ID:%d @%.2fs):\t%s\n',obj.name,obj.eventID,obj.time,obj.info);
        end
    end
end

repo: OpenMAS