text:
% test of Laplacian pyramid

n = 256;
name = 'lena';
M = load_image(name,n);
Jmin = 3;
MW = perform_pyramid_transform_simoncelli(M, Jmin);

save_image = 0;

for i=1:4
    subplot(2,2,i);
    imagesc(MW{i});
    axis image;
    axis off;
    if save_image
        imwrite(rescale(MW{i}), sprintf('MW%d.png', i), 'png');
    end
end
colormap gray(256);

if save_image
    imwrite(rescale(M), sprintf('M.png', i), 'png');
end

author: gpeyre

repo: matlab-toolboxes


text:
function r8vec_convolution_circ_test ( )

%*****************************************************************************80
%
%% R8VEC_CONVOLUTION_CIRC_TEST tests R8VEC_CONVOLUTION_CIRC
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 4;

  x = [ 1.0, 2.0, 3.0, 4.0 ];
  y = [ 1.0, 2.0, 4.0, 8.0 ];
  z_correct = [ 37.0, 44.0, 43.0, 26.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_CONVOLUTION_CIRC_TEST\n' );
  fprintf ( 1, '  R8VEC_CONVOLUTION_CIRC computes the circular convolution\n' );
  fprintf ( 1, '  of two vectors.\n' );

  r8vec_print ( n, x, '  The factor X:' );
  r8vec_print ( n, y, '  The factor Y:' );

  z = r8vec_convolution_circ ( n, x, y );

  r8vec_print ( n, z, '  The circular convolution z = xCCy:' );

  r8vec_print ( n, z_correct, '  Correct answer:' );

  return
end

author: johannesgerer

repo: jburkardt-m


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

%%
clear; close all; clc;

%% Syntax
% |[varargout]=seedIndex2triangulation(F,V,seedIndex);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function varargout = goertzel(x,fs,ft,N,varargin)
% GOERTZEL   Goertzel Algorithm.
%
%  [RE,IM] = GOERTZEL(X,FS,FT,N) calculates the goertzel algorithm
%  for X, considering a sampling frequency FS, the target
%  frequency FT and using a block of size N.
%
%  [MAG,PHASE] = GOERTZEL(...,'mag') returns the magnitude and
%  phase information of the transform.
%
%  [...] = GOERTZEL(...,'overlap',OL) will perform the calculation
%  on blocks that overlap. The amount of overlapping is defined by
%  OL, and can be set from 0 to 1.
%

[mag,OL] = parse_inputs(varargin{:});

if size(x,1)~=1 & size(x,2)~=1
    error('X must be a vector.')
end

flag_rot = 0;
if size(x,2) == 1
    flag_rot = 1;
    x = x';
end

n = length(x);
k = 1*(0.5+N*ft/fs);
w = 2*pi*k/N;
c = cos(w);
s = sin(w);
coef = 2*c;
ADV = max(round((1-OL)*N),1);

re = zeros(1,floor(n/ADV));
im = re;

for j = 1 : floor((n-N)/ADV)
    q1 = 0;
    q2 = 0;
    for i = 1 : N
        q0 = coef*q1 - q2 + x(i+(j-1)*ADV);
        q2 = q1;
        q1 = q0;
    end
    re(j) = q1 - q2*c;
    im(j) = q2*s;
end

if mag == 0
    varargout{1} = re;
    varargout{2} = im;
else
    ma = sqrt(re.^2 + im.^2);
    ph = angle(re+i*im);

    varargout{1} = ma;
    varargout{2} = ph;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [mag,OL] = parse_inputs(varargin)

mag = 0;
OL = 0;
flag_par = 0;

for i = 1 : nargin
    if flag_par == 1
        flag_par = 0;
        continue
    elseif strcmp(varargin{i},'mag')
        mag = 1;
    elseif strcmp(varargin{i},'overlap')
        if nargin == i
            error('OL not defined.')
        end
        flag_par = 1;
        OL = varargin{i+1};
    else
        error('Unknown parameter.')
    end
end

if OL>1 | OL<0
    error('OL should be set inside the interval [0,1].')
end

author: Sable

repo: mcbench-benchmarks


text:
function [ t, wts ] = cgqfs ( nt, kind, alpha, beta, lo )

%*****************************************************************************80
%
%% CGQFS computes knots and weights of a Gauss quadrature formula.
%
%  Discussion:
%
%    This routine computes the knots and weights of a Gauss quadrature
%    formula with:
%
%    * a classical weight function with default values for a,b
%    * only simple knots
%    * optionally print knots and weights and a check of the moments
%
%    Use routine EIQFS to evaluate a quadrature formula computed by
%    this routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, integer LO, selects the action.
%    > 0, compute and print knots and weights.  Print moments check.
%    = 0, compute knots and weights.
%    < 0, compute and print knots and weights.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Compute the Gauss quadrature formula.
%
  [ t, wts ] = cdgqf ( nt, kind, alpha, beta );
%
%  Exit if no print required.
%
  if ( lo == 0 )
    return
  end

  key = 1;
  mop = 2 * nt;
  m = mop + 1;
  mmex = max ( 2 * nt + 3, 1 );
%
%  All knots have multiplicity = 1.
%
  mlt = zeros(nt,1);
  mlt(1:nt) = 1;
%
%  NDX(I) = I.
%
  ndx = ( 1 : nt );

  w = zeros(m,1);

  chkqfs ( t, wts, mlt, nt, nt, ndx, key, mop, mmex, ...
    kind, alpha, beta, lo, w );

  return
end


text:
function imdb = getImdbDCFNetUn_entropy(varargin)
opts = [];
opts.dataDir = fullfile('..','data');

opts.visualization = true;
opts.output_size = 125;
opts.padding = 2;
opts = vl_argparse(opts, varargin);

% -------------------------------------------------------------------------
%   full dataset:
%           vid2015: 464873   494211
% -------------------------------------------------------------------------
set_name = {'vid2015'};
num_all_frame = 494211;
%% Be careful!!! It takes a HUGE RAM for fast speed!
imdb.images.set = int8(ones(1, num_all_frame));
imdb.images.set(randperm(num_all_frame, 1000)) = int8(2);
imdb.images.images = zeros(opts.output_size, opts.output_size, 3, num_all_frame, 'uint8');
imdb.images.up_index = zeros(1, num_all_frame, 'double'); % The farthest frame can it touch
now_index = 0;

% -------------------------------------------------------------------------
%                                                                   VID2015
% -------------------------------------------------------------------------
if any(strcmpi(set_name, 'vid2015'))
    disp('VID2015 Data:');
    if exist('vid_2015_seg.mat', 'file')
        load('vid_2015_seg.mat');%% use dataPreprocessing;
    else
        error('You should generate <vid_2015_seg.mat> according ''dataPreprocessing'' at first.')
    end
    videos = seg;
    n_videos = numel(videos);
    for  v = 1:n_videos
        video = videos{v};fprintf('%3d / %3d \n', v, n_videos);

        img_files = video.path;
        img_first = vl_imreadjpeg(img_files(1));
        im_first = uint8(img_first{1});
        [H, W, ~] = size(im_first);
        img_num = length(img_files);
        % box:(x1, y1, x2, y2)
        % bbox position selection
        boxes = zeros(25,4);
        entropy = zeros(25,1);
        boxes(1,:) = round([3*W/12, 3*H/12, 5*W/12, 5*H/12]);
        boxes(2,:) = round([4*W/12, 3*H/12, 6*W/12, 5*H/12]);
        boxes(3,:) = round([5*W/12, 3*H/12, 7*W/12, 5*H/12]);
        boxes(4,:) = round([6*W/12, 3*H/12, 8*W/12, 5*H/12]);
        boxes(5,:) = round([7*W/12, 3*H/12, 9*W/12, 5*H/12]);

        boxes(6,:) = round([3*W/12, 4*H/12, 5*W/12, 6*H/12]);
        boxes(7,:) = round([4*W/12, 4*H/12, 6*W/12, 6*H/12]);
        boxes(8,:) = round([5*W/12, 4*H/12, 7*W/12, 6*H/12]);
        boxes(9,:) = round([6*W/12, 4*H/12, 8*W/12, 6*H/12]);
        boxes(10,:) = round([7*W/12, 4*H/12, 9*W/12, 6*H/12]);

        boxes(11,:) = round([3*W/12, 5*H/12, 5*W/12, 7*H/12]);
        boxes(12,:) = round([4*W/12, 5*H/12, 6*W/12, 7*H/12]);
        boxes(13,:) = round([5*W/12, 5*H/12, 7*W/12, 7*H/12]);
        boxes(14,:) = round([6*W/12, 5*H/12, 8*W/12, 7*H/12]);
        boxes(15,:) = round([7*W/12, 5*H/12, 9*W/12, 7*H/12]);

        boxes(16,:) = round([3*W/12, 6*H/12, 5*W/12, 8*H/12]);
        boxes(17,:) = round([4*W/12, 6*H/12, 6*W/12, 8*H/12]);
        boxes(18,:) = round([5*W/12, 6*H/12, 7*W/12, 8*H/12]);
        boxes(19,:) = round([6*W/12, 6*H/12, 8*W/12, 8*H/12]);
        boxes(20,:) = round([7*W/12, 6*H/12, 9*W/12, 8*H/12]);

        boxes(21,:) = round([3*W/12, 7*H/12, 5*W/12, 9*H/12]);
        boxes(22,:) = round([4*W/12, 7*H/12, 6*W/12, 9*H/12]);
        boxes(23,:) = round([5*W/12, 7*H/12, 7*W/12, 9*H/12]);
        boxes(24,:) = round([6*W/12, 7*H/12, 8*W/12, 9*H/12]);
        boxes(25,:) = round([7*W/12, 7*H/12, 9*W/12, 9*H/12]);
        for i = 1:25
            entropy(i) = compute_entropy(im_first, boxes(i,:));
        end
        [~, id] = sort(entropy, 'descend');
        % box = round([5*W/12, 5*H/12, 7*W/12, 7*H/12]);
        box = boxes(id(1),:);

        target_pos = [box(1)+box(3), box(2)+box(4)]/2;
        target_sz = [box(3)-box(1), box(4)-box(2)];
        tracklet = KCFtracker(img_files, target_pos, target_sz);
        bbox = [tracklet(:,1)-target_sz(1)/2, tracklet(:,2)-target_sz(2)/2, tracklet(:,1)+target_sz(1)/2, tracklet(:,2)+target_sz(2)/2];

        % bbox = repmat(box, [img_num, 1]);

        im_bank = vl_imreadjpeg(img_files, 'Pack', 'Resize', [H, W], 'numThreads', 32);
        n_frames = size(bbox,1);
        imdb.images.images(:,:,:,now_index+(1:n_frames)) = uint8(...
            imcrop_pad(im_bank{1}, bbox, opts.padding, opts.output_size([1,1])));
        imdb.images.up_index(now_index+(1:n_frames)) = (n_frames:-1:1)-1;
        imdb.images.set(now_index+n_frames) = 4; %should not be selected as x.
        now_index = now_index + n_frames;
    end %%end v
end %%end VID2017

% dataMean = single(mean(imdb.images.images,4));
dataMean(1,1,1:3) = single([123,117,104]);
imdb.images.data_mean(1, 1, 1:3) = dataMean;
imdb.meta.sets = {'train', 'val'} ;

end %%end function

function patch_entropy = compute_entropy(image, box)

% box: [x1, y1, x2, y2]
pos = [box(1)+box(3), box(2)+box(4)]/2;
sz = ([box(3)-box(1), box(4)-box(2)])*(1); % consider small padding
xs = round(pos(1) + (1:sz(1)) - sz(1)/2);
ys = round(pos(2) + (1:sz(2)) - sz(2)/2);
%extract image
im_crop = image(ys, xs, :);
patch_entropy = entropy(im_crop);

end

repo: UDT


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

dy = 0.15;
dx = 0.15 ;
dz = 0.30 ;
ni = 150;

%a.Longitude = (a.Longitude - mean(a.Longitude))*111;
%a.Latitude = (a.Latitude - mean(a.Latitude))*111;
xc = 0; yc = 0; zc =0 ;
x0 = min(a.Longitude); x1 = max(a.Longitude);
y0 = min(a.Latitude); y1 = max(a.Latitude);
z0 = min(a.Depth); z1 = max(a.Depth);
xv = x0:dx:x1;yv= y0:dy:y1;zv = z0:dz:z1;
bvg = ones(length(x0:dx:x1),length(y0:dy:y1),length(z0:dz:z1));
ra = ones(length(x0:dx:x1),length(y0:dy:y1),length(z0:dz:z1));
itotal = length(x0:dx:x1)*length(y0:dy:y1)*length(z0:dz:z1)

wai = waitbar(0,' Please Wait ...  ');
set(wai,'NumberTitle','off','Name','b-value grid - percent done');;
drawnow

allcount = 0;
i2 = 0;
for x = x0:dx:x1
    xc = xc+1;
    for y = y0:dy:y1
        yc = yc+1;
        for z = z0:dz:z1
            zc = zc+1;
            allcount = allcount + 1.;
            i2 = i2+1;

            % calculate distance from center point and sort wrt distance
            l = sqrt(((a.Longitude - x)).^2 + ((a.Latitude - y).^2 + ((a.Depth - z)).^2)) ;
            [s,is] = sort(l);
            b = a(is(:,1),:) ;       % re-orders matrix to agree row-wise

            % take first ni points
            b = b(1:ni,:);      % new data per grid point (b) is sorted in distance

            % call the b-value function
            [bv magco stan av me mer me2,  pr] =  bvalca3(b,inb1,inb2);
            l2 = sort(l);
            b2 = b;
            %if inb2 ==  1
            %l = b(:,6) >= magco;
            %b2 = b(l,:);
            %end
            %[av2 bv2 stan2 ] =  bmemag(b2);
            bvg(xc,yc,zc) = bv;
            ra(xc,yc,zc) = l2(ni);
            waitbar(allcount/itotal)
        end   % for z
        zc = 0;
    end   % for y
    yc = 0;
end   % for x
xc = 0;

close(wai)

bv2 = bvg;
bvg = bv2;
l = ra > 2.500;
bvg(l)=nan;
figure
[X,Y,Z] = meshgrid(yv,xv,zv);
zs = [z0:4*dz:z1]
sl = slice(X,Y,Z,bvg,[x0+3*dx],[y0+4*dy],zs)
clf
sl1 = slice(X,Y,Z,bvg,x1-6*dx,y1-6*dy,[ z0+4*dz 9 14])
hold on
%plot3(a.Latitude,a.Longitude,a.Depth,'k.','MarkerSize',2)
rotate3d on
caxis([1.8 2.8])
set(gca,'XLim',[-3 3 ],'xgrid','off')
set(gca,'YLim',[-3 3 ],'ygrid','off')
set(gca,'ZLim',[  -1 20 ],'zgrid','off')
colormap(h)
shading interp
cob = colorbar('vert')
set(cob,'TickDir','out','pos',[0.8 0.3 0.07 0.3])
set(gca,'Box','on','vis','on')
tmp = ra*nan;
tmp(1,1,1) = 0;
tmp(1,1,2) = 1;
hold on
sl = slice(X,Y,Z,tmp,x1-6*dx,y1-6*dy,[ z0+4*dz 9 14 ])
set(sl(:),'EdgeColor','w')
caxis([1.8 2.8])
view([-36 10])

author: CelsoReyes


text:
function val = get_params(CPD, name)
% GET_PARAMS Get the parameters (fields) for a tabular_CPD object
% val = get_params(CPD, name)
%
% The following fields can be accessed
%
% cpt       - the CPT
%
% e.g., CPT = get_params(CPD, 'cpt')

switch name
 case 'cpt',      val = CPD.CPT;
 case 'tree',     val = CPD.tree;
 otherwise,
  error(['invalid argument name ' name]);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%{

This code is to render a Mesh given a 3x4 camera matrix with an image resolution widthxheight. The rendering result is an ID map for facets, edges and vertices. This can usually used for occlusion testing in texture mapping a model from an image, such as the texture mapping in the following two papers.

--Jianxiong Xiao http://mit.edu/jxiao/

Citation:

[1] J. Xiao, T. Fang, P. Zhao, M. Lhuillier, and L. Quan
Image-based Street-side City Modeling
ACM Transaction on Graphics (TOG), Volume 28, Number 5
Proceedings of ACM SIGGRAPH Asia 2009

[2] J. Xiao, T. Fang, P. Tan, P. Zhao, E. Ofek, and L. Quan
Image-based Facade Modeling
ACM Transaction on Graphics (TOG), Volume 27, Number 5
Proceedings of ACM SIGGRAPH Asia 2008

%}

clear
clc
close all

compile
load('data.mat')

result = RenderMex(P, img_width, img_height, vertex, edge, face)';

close all
imagesc(result)
axis equal
axis tight
max(max(result))

author: jianxiongxiao

repo: ProfXkit


text:
function output = picker(condition, input1, input2)
    if condition
        output = input1;
    else
        output = input2;
    end
end

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
addpath('tools')
more off
close all
clear all

% Load laser scans and robot poses.
load("../data/laser")
%laser = read_robotlaser('../data/csail.log');

% Extract robot poses: Nx3 matrix where each row is in the form: [x y theta]
poses = [laser.pose];
poses = reshape(poses,3,size(poses,2)/3)';

% Initial cell occupancy probability.
prior = 0.50;
% Probabilities related to the laser range finder sensor model.
probOcc = 0.9;
probFree = 0.35;

% Map grid size in meters. Decrease for better resolution.
gridSize = 0.1;

% Set up map boundaries and initialize map.
border = 30;
robXMin = min(poses(:,1));
robXMax = max(poses(:,1));
robYMin = min(poses(:,2));
robYMax = max(poses(:,2));
mapBox = [robXMin-border robXMax+border robYMin-border robYMax+border];
offsetX = mapBox(1);
offsetY = mapBox(3);
mapSizeMeters = [mapBox(2)-offsetX mapBox(4)-offsetY];
mapSize = ceil([mapSizeMeters/gridSize]);

% Used when updating the map. Assumes that prob_to_log_odds.m
% has been implemented correctly.
logOddsPrior = prob_to_log_odds(prior);

% The occupancy value of each cell in the map is initialized with the prior.
map = logOddsPrior*ones(mapSize);
disp('Map initialized. Map size:'), disp(size(map))

% Map offset used when converting from world to map coordinates.
offset = [offsetX; offsetY];

% Main loop for updating map cells.
% You can also take every other point when debugging to speed up the loop (t=1:2:size(poses,1))
for(t=1:size(poses,1))
%for(t=1:50)
	t
	% Robot pose at time t.
	robPose = [poses(t,1);poses(t,2);poses(t,3)];

	% Laser scan made at time t.
	sc = laser(1,t);
	% Compute the mapUpdate, which contains the log odds values to add to the map.
	[mapUpdate, robPoseMapFrame, laserEndPntsMapFrame] = inv_sensor_model(map, sc, robPose, gridSize, offset, probOcc, probFree);

	mapUpdate -= logOddsPrior*ones(size(map));
	% Update the occupancy values of the affected cells.
	map += mapUpdate;

	% Plot current map and robot trajectory so far.
        plot_map(map, mapBox, robPoseMapFrame, poses, laserEndPntsMapFrame, gridSize, offset, t);
endfor

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function [c1,c2]= singlePtX(p1,p2,bounds,Ops)
% function [c1,c2] = singlePtXover(p1,p2,bounds,Ops)
% p1      - the first parent ( [solution string function value] )
% p2      - the second parent ( [solution string function value] )
% bounds  - the bounds matrix for the solution space
% Ops     - Options matrix for simple crossover [gen #SimpXovers].

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% 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 1, 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. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
sz=size(p1,2)-1;
cut = round(rand*(sz-1)+0.5); %Generate random cut point U(1,n-1)
pm1=p1(1:sz);
pm2=p2(1:sz);
c1=p1;
c2=p2;
c1(1:cut)=p1(1:cut);
c2(1:cut)=p2(1:cut);
for i=1:cut
  pm1=strrep(pm1,p2(i),-1);
  pm2=strrep(pm2,p1(i),-1);
end
c1((cut+1):sz)=p2(find(pm2>0));
c2((cut+1):sz)=p1(find(pm1>0));

author: HuangCongQing


text:
function [C, dA] = mcml_grad(x, X, P)
%MCML_GRAD Computes MCML cost function and gradient
%
%   [C, dA] = mcml_grad(x, X, P)
%
% Computes MCML cost function and gradient.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    % Decode solution
    [n, d] = size(X);
    A = reshape(x, [d d]);

    % Compute conditional probabilities for current solution
    D = zeros(n, n);
    for i=1:n
        diffX = bsxfun(@minus, X(i,:), X(i + 1:end,:));
        dist = sum((diffX * A) .* diffX, 2);
        D(i + 1:end, i) = dist;
        D(i, i + 1:end) = dist';
    end
    Q = exp(-D);
    Q(1:n+1:end) = 0;
    Q = bsxfun(@rdivide, Q, sum(Q, 2));
    Q = max(Q, realmin);

    % Compute cost function
    C = sum(P(:) .* log(P(:) ./ Q(:)));

    % Compute gradient with respect to A
    if nargin > 1
        dA = zeros(d, d);
        PQ = P - Q;
        for i=1:n
            diffX = bsxfun(@minus, X(i,:), X);
            dA = dA + bsxfun(@times, PQ(i,:), diffX') * diffX;
        end
    end
    dA = dA(:);

repo: cluster


text:
function [x, F, exitflag] = NewtonRaphson_octave(fun, x0, options)
% NEWTONRAPHSON Solve set o non-linear equations using Newton-Raphson method.
% this version works in Octave

% Copyright (C) 2021 Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% This code is modified from the original found at:
% https://github.com/mikofski/NewtonRaphson

% Copyright (c) 2021, Mark Mikofski
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice, this
%    list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
%    this list of conditions and the following disclaimer in the documentation
%    and/or other materials provided with the distribution.
%
% 3. Neither the name of the copyright holder nor the names of its
%    contributors may be used to endorse or promote products derived from
%    this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% [X, RESNORM, F, EXITFLAG] = NEWTONRAPHSON(FUN, X0, OPTIONS)
% FUN is a function handle that returns a vector of residuals equations, F,
% and takes a vector, x, as its only argument. When the equations are
% solved by x, then F(x) == zeros(size(F(:), 1)).
%
% Optionally FUN may return the Jacobian, Jij = dFi/dxj, as an additional
% output. The Jacobian must have the same number of rows as F and the same
% number of columns as x. The columns of the Jacobians correspond to d/dxj and
% the rows correspond to dFi/d.
%
%   EG:  J23 = dF2/dx3 is the 2nd row ad 3rd column.
%
% If FUN only returns one output, then J is estimated using a center
% difference approximation,
%
%   Jij = dFi/dxj = (Fi(xj + dx) - Fi(xj - dx))/2/dx.
%
% NOTE: If the Jacobian is not square the system is either over or under
% constrained.
%
% X0 is a vector of initial guesses.
%
% OPTIONS is a structure of solver options created using OPTIMSET.
% EG: options = optimset('TolX', 0.001).
%
% The following options can be set:
% * OPTIONS.TOLFUN is the maximum tolerance of the norm of the residuals.
%   [1e-6]
% * OPTIONS.TOLX is the minimum tolerance of the relative maximum stepsize.
%   [1e-6]
% * OPTIONS.MAXITER is the maximum number of iterations before giving up.
%   [100]
%
% X is the solution that solves the set of equations within the given tolerance.
% RESNORM is norm(F) and F is F(X). EXITFLAG is an integer that corresponds to
% the output conditions, OUTPUT is a structure containing the number of
% iterations, the final stepsize and exitflag message and JACOB is the J(X).
%
% See also OPTIMSET, OPTIMGET, FMINSEARCH, FZERO, FMINBND, FSOLVE, LSQNONLIN
%
% References:
% * http://en.wikipedia.org/wiki/Newton's_method
% * http://en.wikipedia.org/wiki/Newton's_method_in_optimization
% * 9.7 Globally Convergent Methods for Nonlinear Systems of Equations 383,
%   Numerical Recipes in C, Second Edition (1992),
%   http://www.nrbook.com/a/bookcpdf.php
% Version 0.5
% * allow sparse matrices, replace cond() with condest()
% * check if Jstar has NaN or Inf, return NaN or Inf for cond() and return
%   exitflag: -1, matrix is singular.
% * fix bug: max iteration detection and exitflag reporting typos
% Version 0.4
% * allow lsq curve fitting type problems, IE non-square matrices
% * exit if J is singular or if dx is NaN or Inf
% Version 0.3
% * Display RCOND each step.
% * Replace nargout checking in funwrapper with ducktypin.
% * Remove Ftyp and F scaling b/c F(typx)->0 & F/Ftyp->Inf!
% * User Numerical Recipies minimum Newton step, backtracking line search
%   with alpha = 1e-4, min_lambda = 0.1 and max_lambda = 0.5.
% * Output messages, exitflag and min relative step.
% Version 0.2
% * Remove `options.FinDiffRelStep` and `options.TypicalX` since not in MATLAB.
% * Set `dx = eps^(1/3)` in `jacobian` function.
% * Remove `options` argument from `funwrapper` & `jacobian` functions
%   since no longer needed.
% * Set typx = x0; typx(x0==0) = 1; % use initial guess as typx, if 0 use 1.
% * Replace `feval` with `evalf` since `feval` is builtin.
%% initialize
% There are no argument checks!
x0 = x0(:); % needs to be a column vector
% set default options

defaultopt = struct('TolX', 1e-12, 'TolFun', 1e-6, 'MaxIter', 1000);%, 'Display', 'off');

%% get options
TOLX = optimget(options, 'TolX', defaultopt.TolX);
TOLFUN = optimget(options, 'TolFun', defaultopt.TolFun);
MAXITER = optimget(options, 'MaxIter', defaultopt.MaxIter);

%TYPX = max(abs(x0), 1); % x scaling value, remove zeros
ALPHA = 1e-4; % criteria for decrease
MIN_LAMBDA = 0.1; % min lambda
MAX_LAMBDA = 0.5; % max lambda
%% set scaling values
% TODO: let user set weights
%weight = ones(numel(fun(x0)),1);
%J0 = weight*(1./TYPX'); % Jacobian scaling matrix
%% check initial guess
x = x0; % initial guess
F = fun(x); % evaluate initial guess
nf = length(F);
J = jacobian(fun, x, nf, F);
%Jstar = J./J0; % scale Jacobian
%if any(isnan(Jstar(:))) || any(isinf(Jstar(:)))
if any(isnan(J(:))) || any(isinf(J(:)))
    exitflag = -1; % matrix may be singular
else
    exitflag = 1; % normal exit
end
resnorm = norm(F, Inf); % calculate norm of the residuals
resnorm0 = 100*resnorm;dx = zeros(size(x0)); % dummy values
%% solver
Niter = 0; % start counter
lambda = 1; % backtracking
while (resnorm>TOLFUN || lambda<1) && exitflag>=0 && Niter<=MAXITER
    if lambda==1
        %% Newton-Raphson solver
        Niter = Niter+1; % increment counter
        % update Jacobian, only if necessary
        if resnorm/resnorm0 > .2
            J = jacobian(fun, x, nf, F);
%            Jstar = J./J0; % scale Jacobian
%            if any(isnan(Jstar(:))) || any(isinf(Jstar(:)))
            if any(isnan(J(:))) || any(isinf(J(:)))
                exitflag = -1; % matrix may be singular
                break
            end
        end
%        dx_star = -Jstar\F; % calculate Newton step
%        % NOTE: use isnan(f) || isinf(f) instead of STPMAX
%        dx = dx_star.*TYPX; % rescale x
        if rcond(J) <= eps; dx = pinv(J) * -F; else dx = -J\F; end
        g = F'*J;%star; % gradient of resnorm
        slope = g*dx;%_star; % slope of gradient
        fold = F'*F; % objective function
        xold = x; % initial value
        lambda_min = TOLX/max(abs(dx)./max(abs(xold), 1));
    end
    if lambda<lambda_min
        exitflag = 2; % x is too close to XOLD
        break
    elseif any(isnan(dx)) || any(isinf(dx))
        exitflag = -1; % matrix may be singular
        break
    end
    x = xold+dx*lambda; % next guess
    F = fun(x); % evaluate this guess
    f = F'*F; % new objective function
    %% check for convergence
    lambda1 = lambda; % save previous lambda
    if f>fold+ALPHA*lambda*slope
        if lambda==1
            lambda = -slope/2/(f-fold-slope); % calculate lambda
        else
            A = 1/(lambda1 - lambda2);
            B = [1/lambda1^2,-1/lambda2^2;-lambda2/lambda1^2,lambda1/lambda2^2];
            C = [f-fold-lambda1*slope;f2-fold-lambda2*slope];
            coeff = num2cell(A*B*C);
            [a,b] = coeff{:};
            if a==0
                lambda = -slope/2/b;
            else
                discriminant = b^2 - 3*a*slope;
                if discriminant<0
                    lambda = MAX_LAMBDA*lambda1;
                elseif b<=0
                    lambda = (-b+sqrt(discriminant))/3/a;
                else
                    lambda = -slope/(b+sqrt(discriminant));
                end
            end
            lambda = min(lambda,MAX_LAMBDA*lambda1); % minimum step length
        end
    elseif isnan(f) || isinf(f)
        % limit undefined evaluation or overflow
        lambda = MAX_LAMBDA*lambda1;
    else
        lambda = 1; % fraction of Newton step
    end
    if lambda<1
        lambda2 = lambda1;f2 = f; % save 2nd most previous value
        lambda = max(lambda,MIN_LAMBDA*lambda1); % minimum step length
        continue
    end
    resnorm0 = resnorm; % old resnorm
    resnorm = norm(F, Inf); % calculate new resnorm
end
%% output
% output.iterations = Niter; % final number of iterations
% output.stepsize = dx; % final stepsize
% output.lambda = lambda; % final lambda
if Niter>=MAXITER
    exitflag = 0;
%     output.message = 'Number of iterations exceeded OPTIONS.MAXITER.';
% elseif exitflag==2
%     output.message = 'May have converged, but X is too close to XOLD.';
% elseif exitflag==-1
%     output.message = 'Matrix may be singular. Step was NaN or Inf.';
% else
%     output.message = 'Normal exit.';
end
% jacob = J;
end
function J = jacobian(fun, x, nf, funx)
% estimate J
dx = eps^(1/3); % finite difference delta
nx = numel(x); % degrees of freedom
if nargin <4
    funx = fun(x);
    if nargin <3
        nf = numel(funx);
    end % number of functions
end
J = zeros(nf,nx); % matrix of zeros
for n = 1:nx
    % create a vector of deltas, change delta_n by dx
    delta = J(:, n); delta(n) = dx;
    dF = fun(x+delta)-funx;%-delta); % delta F
    J(:, n) = dF(:)/dx;%/2; % derivatives dF/d_n
end
end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [returnVar,msg] = ReRefFilFile(fname,nbChan,refChan,rerefChan)

% USAGE:
%     ReRefFilFile(fbasename, nbChan, refChan, rerefChan)
%     This function substracts one channel or the median of a series of
%     channels (the 'reference' channel) from other channels.
%
%     Note: this uses 1-based indexing, not zero-based, so you have to take
%     the channel number in neuroscope and add one to it, e.g. channel 72
%     in neuroscope is channel 73 here.
%
% INPUTS:
%     fname: fil file name. Could be also a file basename not finishing
%     with '.fil' and the program will be then executed for all the fil
%     files begining with this file basename.
%     nbChan: total number of channels in fil file
%     refChan: reference channel (could be a vector, in this case the median will be the new reference)
%     rerefChan: vector of channels to rereference
%
% Adrien Peyrache 2013

if ~strcmp(fname(end-2:end),'fil')
    datFiles = dir([fname '*.fil']);
    for ii=1:length(datFiles)
        ReRefFilFile(datFiles(ii).name,nbChan,refChan,rerefChan)
    end

else
    fprintf('ReReferencing %s\n',fname)
    try
        infoFile = dir(fname);

        chunk = 1e6;
        nbChunks = floor(infoFile.bytes/(nbChan*chunk*2));
        warning off
        if nbChunks==0
            chunk = infoFile.bytes/(nbChan*2);
        end

        for ix=0:nbChunks-1
            m = memmapfile(fname,'Format','int16','Offset',ix*chunk*nbChan*2,'Repeat',chunk*nbChan,'writable',true);
            d = m.Data;
            d = reshape(d,[nbChan chunk]);
            ref = d(refChan,:);
            if length(refChan)>1
                ref = median(double(ref));
            end
            ref = repmat(ref,[length(rerefChan) 1]);
            d(rerefChan,:) = d(rerefChan,:)-int16(ref);

            m.Data = d(:);
            clear d m
        end
        %close(h)

        newchunk = infoFile.bytes/(2*nbChan)-nbChunks*chunk;

        if newchunk
            m = memmapfile(fname,'Format','int16','Offset',nbChunks*chunk*nbChan*2,'Repeat',newchunk*nbChan,'writable',true);
            d = m.Data;
            d = reshape(d,[nbChan newchunk]);
            ref = d(refChan,:);
            if length(refChan)>1
                ref = median(ref);
            end
            ref = repmat(ref,[length(rerefChan) 1]);
            d(rerefChan,:) = d(rerefChan,:)-int16(ref);
            m.Data = d(:);
         clear d m
        end
        warning on
        returnVar = 1;
        msg = '';

    catch
        fprintf(['Error occurred in processing ' fname '. File not rereferenced.\n']);
        keyboard
        returnVar = 0;
        msg = lasterr;
    end
    clear m
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ mcmap, cdata_mapped ] = multicmap( h, cmaps, clims )
% MULTICMAP Apply multiple colormaps to image objects
%
%   [MCMAP,CDATA]=MULTICMAP(H,CMAPS,CLIMS) given handles to
%   image objects plotted in the current axes, along with
%   specification of individual colormaps for the image
%   objects, this function combines the individual colormaps
%   into a single multi-colormap matrix, and computes re-mapped
%   CData matrices for individual image objects, so that they
%   point to their respective colormaps within the multi-colormap.
%
%   This function assumes that CData of each image object
%   is a 2D matrix of indexes to the figures colormap
%   (which is the default behavior in MATLAB), as opposed
%   to CData being a 3D matrix of true color RGB image values.
%
%   Not that this routine is a generalization of examples
%   presented in [1].
%
%   Known limitations: currently multiple subplots in a single
%                      figure are not supported.
%
%   Inputs
%           H numeric array (or structure) containing handles
%             to image objects plotted in the current axes
%
%           CMAPS column cell array (or structure) of colormaps
%                 for the image objects specified in H. Each
%                 colormap should conform to MATLAB's colormap
%                 specifications, i.e., it should be an m-by-3
%                 matrix of real numbers between 0.0 and 1.0.
%                 Each row of a given colormap matrix is an RGB
%                 vector that defines one color. See "doc colormap"
%                 for further information regarding MATLAB's colormaps.
%
%           CLIMS m-by-2 matrix, or structure of two element
%                 row vectors, each containing lower and upper
%                 limits for each image object's CData indexes.
%                 See "doc imagesc" for detailed diagram and
%                 explanation of CLIMS.
%
%   Outputs
%           MCMAP m-by-3 colormap matrix composed of individual
%                 colormaps supplied in CMAPS.
%
%           CDATA re-mapped CData matrices as members of a cell array
%                 for individual image objects specified in H.
%
%   References
%           [1] MathWorks: Product Support,
%               "1215 - Using Multiple Colormaps in a Single Figure"
%               url: http://www.mathworks.com/support/tech-notes/1200/1215.html
%
%   Example
%            clear all; close all; clc;
%
%            % Read the sample sample images
%            data.earth = load( 'earth' );
%            data.penny = load( 'penny' );
%
%            % Get indexed images into variables
%            photo.earth = data.earth.X;
%            photo.penny = data.penny.P;
%
%            % Get dimensions of each image
%            dims.earth = size( photo.earth );
%            dims.penny = size( photo.penny );
%
%            % Define x- and y-points for each image
%            x.earth = [ 1:dims.earth(1) ];
%            x.penny = [ 1:dims.penny(1) ] + round((dims.earth(1)-dims.penny(1))/2);
%            y.earth = 1:dims.earth(2);
%            y.penny = [ 1:dims.penny(2) ] + round((dims.earth(2)-dims.penny(2))/2);
%
%            % Define colormaps for each image like so:
%            maps.earth = data.earth.map;
%            maps.penny = hot(64);
%
%            % Define transparency for visible parts of each image
%            alpha.earth = 1;        % fully opaque
%            alpha.penny = 0.7;      % 30% transparency
%
%            % Define which parts of the second image are to be fully transparent (invisible)
%            AlphaData.penny = ones( size( photo.penny ) ) * alpha.penny;
%            AlphaData.penny( photo.penny<5*min(photo.penny(:)) ) = 0;
%
%            % Plot images with their respective colormaps
%            hfig = figure( 'Position', [ 400 10 600 600 ], 'PaperPositionMode', 'auto', 'color', 'w' );
%
%            % Make sure both images get retained in the current axes
%            hold on;
%
%            % Plot images and retain handles to the image objects
%            h.earth = image( x.earth, y.earth, photo.earth );
%            h.penny = image( x.penny, y.penny, photo.penny );
%
%            % Apply transparency settings
%            set( h.earth, 'AlphaData', alpha.earth );
%            set( h.penny, 'AlphaData', AlphaData.penny );
%
%            % Apply axes limits
%            xlim( [ min(struct2array(x)) max(struct2array(x)) ] );
%            ylim( [ min(struct2array(y)) max(struct2array(y)) ] );
%
%            % Make sure our images are not up-side-down
%            axis ij square off
%
%            % Apply colormaps to their respective image objects
%            multicmap( h, maps );
%
%   See also EXAMPLE_TWO_IMAGES, EXAMPLE_THREE_IMAGES, EXAMPLE_SPEECH.

%   Author: Kamil Wojcicki, UTD, February 2012.

    % check for correct number of input arguments
    if nargin<2 || nargin>3
        error( sprintf('Incorrect number of input arguments.\nType "help %s" for usage help.\n', mfilename) );
    end

    % if image object handles were passed-in as values for members
    % of a structure, then convert the structure to a numeric array
    if isstruct( h )
        h = struct2array( h );
    end

    % determine number of image objects
    N = length( h );

    % if individual colormaps were passed-in as matrices for members
    % of a structure, then convert the structure to a cell array
    if isstruct( cmaps )
        cmaps = struct2cell( cmaps );
    end

    % combine colormaps for individual image objects into
    % a single multi-colormap (as a m-by-3 matrix)
    mcmap = cell2mat( cmaps );

    % determine if data limits were specified
    if nargin==3

        % if the data limits were specified as 1-by-2 arrays
        % for member of a structure, then convert the structure
        % to a cell array of data limits
        if isstruct( clims )
            clims = struct2cell( clims );
        end

        % if the data limits are as 1-by-2 arrays, members of
        % a cell array, then convert the cell array into m-by-2 matrix
        if iscell( clims )
            clims = cell2mat( clims );
        end

        % determine the dimensionality of the data limits matrix
        [ rows, columns ] = size( clims );

        % for each image object there should be lower
        % and upper limit for the data (i.e., exactly 2)
        assert( columns==2 );

        % if only one array of data limit values was supplied
        % then use it for all image objects by replicating
        % it for each
        if rows==1
            clims = repmat( clims, N, 1 );

        % otherwise, make sure that the number of supplied
        % data limits matched the number of image objects
        else
            assert( rows==N );
        end

    end

    % for each image object compute CData re-mapped
    % to index to the new multi-colormap.
    for n = 1:N

        % get the image data (assumed to be indices
        % into the figure's colormap)
        cdata = get( h(n), 'CData' );

        % if climits were supplied, apply them
        if nargin==3
            cdata( cdata<clims(n,1) ) = clims(n,1);
            cdata( cdata>clims(n,2) ) = clims(n,2);
        end

        % figure-out the maximum, minimum and range
        % of the current image object's CData
        cdata_max = max( cdata(:) );
        cdata_min = min( cdata(:) );
        cdata_range = cdata_max - cdata_min;

        % get the size of the colormap for the
        % current image object
        cmap_size = size( cmaps{n}, 1 );

        % re-map the CData for the current image object to its colormap
        cdata_mapped{n} = min( cmap_size, round((cmap_size-1)*(cdata-cdata_min)/cdata_range)+1 );

        % re-map the CData for the current image object to its
        % corresponding colormap within the multi-colormap
        if n>1, cdata_mapped{n} = cdata_mapped{n} + max( cdata_mapped{n-1}(:) ); end

    end

    % if no outputs have been requested, then update data for
    % for the individual image objects to the re-mapped data,
    % also set data limits (for plotting) and apply the new colormap
    if nargout==0

        % for each image object update CData to re-mapped data
        for n = 1:N
            set( h(n), 'CData', cdata_mapped{n} );
        end

        % set CData axis limits (for plotting)
        caxis( [ min(cdata_mapped{1}(:)) max(cdata_mapped{end}(:)) ] );

        % apply the new colormap
        colormap( mcmap );

    end

% EOF

author: optimaltransport


text:
function [thr, h1] = badacostCalibrateCascade(X0, X1, clf, use_trees)
% Classify all the data with the classifier and check the
% costs traces. We compare positive class trace (label 1) with
% the min cost of the positive classes. We use the formula:
%
%  traces(i,:) = -(min_pos_costs - costs(1,:));
%
%  Where traces(i, :) are the traces for the i-th badacost weak learner,
%  costs(1,:) are the costs associated with the negative class at
%  i-th badacost weak-learner and min_pos_cost is the min cost at the
%  i-th weak-learner of the positive classes.

if nargin < 4
  use_trees = 1;
end

% 1. Compute costs traces ...
if use_trees
  nWeaks = size(clf.fids, 2);
else
  nWeaks = length(clf.WEAK_LEARNERS);
end
n0 = size(X0, 1);
n1 = size(X1, 1);
%n  = n0 + n1;
margin_vec0 = zeros(clf.num_classes, n0);
margin_vec1 = zeros(clf.num_classes, n1);
traces0     = zeros(nWeaks, n0);
traces1     = zeros(nWeaks, n1);

for i=1:nWeaks
  if use_trees
    % z is a row vector with the labels
    ids0 = forestInds(X0,clf.thrs(:,i),clf.fids(:,i),clf.child(:,i), 1e5);
    for j=1:n0
      z = clf.hs(ids0(j),i);
      margin_vec0(:,j) = margin_vec0(:,j) + ...
                        (clf.wl_weights(i).*clf.Y(:, z));
    end
  else % General case with other weak learners
    % z is a row vector with the labels
    z  = clf.classify_weak_learner(clf.WEAK_LEARNERS{i}, X0');
    for j=1:n0
      margin_vec0(:,j) = margin_vec0(:,j) + (clf.WEIGHTS(i).*clf.Y(:, z(j)));
    end
  end
  % WARNING: Change to accomodate with theory (2016/11)
%  costs = clf.Cprime' * margin_vec0;
  costs = clf.Cprime * margin_vec0;
  [min_pos_costs , ~] = min(costs(2:end,:));
  traces0(i,:) = -(min_pos_costs - costs(1,:));

  if use_trees
    ids1 = forestInds(X1,clf.thrs(:,i),clf.fids(:,i),clf.child(:,i), 1e5);
    for j=1:n1
      z = clf.hs(ids1(j),i);
      margin_vec1(:,j) = margin_vec1(:,j) + ...
                        (clf.wl_weights(i).*clf.Y(:, z));
    end
  else % general case with other weak learners.
    % z is a row vector with the labels
    z  = clf.classify_weak_learner(clf.WEAK_LEARNERS{i}, X1');
    for j=1:n1
      margin_vec1(:,j) = margin_vec1(:,j) + (clf.WEIGHTS(i).*clf.Y(:, z(j)));
    end
  end
  % WARNING: Change to accomodate with theory (2016/11)
%  costs = clf.Cprime' * margin_vec1;
  costs = clf.Cprime * margin_vec1;
  [min_pos_costs , ~] = min(costs(2:end,:));
  traces1(i,:) = -(min_pos_costs - costs(1,:));
end

% 2. plot the min and max trace at last weak learner of positives and
%    negatives data
[~ , min_neg_trace_index] = min(traces0(nWeaks,:));
[~ , max_neg_trace_index] = max(traces0(nWeaks,:));

[~ , min_pos_trace_index] = min(traces1(nWeaks,:));
[~ , max_pos_trace_index] = max(traces1(nWeaks,:));

h1 = figure;
plot(traces1(:,min_pos_trace_index), 'g-', 'LineWidth', 1.5);
hold on;
plot(traces1(:,max_pos_trace_index), 'g-', 'LineWidth', 1.5);
plot(traces0(:,min_neg_trace_index), 'r-', 'LineWidth', 1.5);
plot(traces0(:,max_neg_trace_index), 'r-', 'LineWidth', 1.5);

% 3. Compute the optimal threshold (for same success rate but with
% less average weak-learners evaluated per feature vector). For doing
% so, we find the trace of the last positive example that went over
% 0 at the last badacost weak-learner.
indices       = find(traces1(nWeaks,:) > 0);
[~, index] = min(traces1(nWeaks, indices));
index = indices(index);

figure(h1);
hold on;
plot(traces1(:,index), 'b-', 'LineWidth', 2);

min_last_traces1_pos = min(traces1(:,index));
thr = min_last_traces1_pos;

hold on;
plot([1 nWeaks], [min_last_traces1_pos min_last_traces1_pos], 'k-', 'LineWidth', 2);
plot([1 nWeaks], [thr thr], 'm-', 'LineWidth', 2);
hold off;

author: MengyangPu

repo: EDTER

save_path: github-repos/MATLAB/MengyangPu-EDTER


text:
function rs = gmi(s, D, eps, NNR, len, Nref)

%tstoolbox/@signal/gmi
%   Syntax:
%     * gmi(s, D, eps, NNR, len, Nref)
%
%   Input arguments:
%     * D -
%     * eps -
%     * NNR -
%     * len -
%     * Nref -
%
%   Generalized mutual information function for a scalar time series
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt
narginchk(6,6);

if (ndim(s) > 1) | (~isreal(data(s)))
	help(mfilename)
	return
end

L = dlens(s, 1);
[y,i] = sort(data(s));
ra(i) =  0:(1/L):(1-1/L);		% create rank values

e = embed(signal(ra'), D, 1);		%'

c = core(gmi(data(e), ceil(rand(Nref,1)*(L-len-D-3))+2, eps, NNR, len));

rs = signal(c, s);      % special constructor calling syntax for working routines
a = getaxis(s, 1);
dl = delta(a);
a = setfirst(a, 0);
rs = setaxis(rs, 1, a);
rs = setyunit(rs, unit);                % gmi values are scalars without unit
rs = addhistory(rs, 'Generalized mutual information function');
rs = addcommandlines(rs, 's = gmi(s', D, eps, NNR, len, Nref);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function LoadFishDataWithoutTS(hfig,i_fish,hdf5_dir,mat_dir)
disp(['loading fish #' num2str(i_fish) '...']);

M_fish_set = getappdata(hfig,'M_fish_set');
fishset = M_fish_set(i_fish);
setappdata(hfig,'fishset',fishset);

%% load data from file

disp(['load fish ' num2str(i_fish) '...']);

if ~exist('hdf5_dir','var') || ~exist('mat_dir','var'),
    data_masterdir = getappdata(hfig,'data_masterdir');
    data_dir = fullfile(data_masterdir,['subject_' num2str(i_fish)]);
    hdf5_dir = fullfile(data_dir,'TimeSeries.h5');
    mat_dir = fullfile(data_dir,'data_full.mat');
end

tic

% load 'data'
load(mat_dir,'data'); % struct with many fields
names = fieldnames(data); % cell of strings
for i = 1:length(names),
    setappdata(hfig,names{i},eval(['data.',names{i}]));
end

% load time series (hdf5 file)
absIX = h5read(hdf5_dir,'/absIX');

setappdata(hfig,'absIX',absIX);

toc

% setappdata(hfig,'isfullfish',1);
end

author: xiuyechen

repo: FishExplorer


text:
clc;
clear all;
close all;
tic;

negfiles = getAllFiles('data\txt_sentoken\neg\');
posfiles = getAllFiles('data\txt_sentoken\pos\');
labels = [-ones(size(negfiles,1),1); ones(size(posfiles,1),1)];

allfiles = [negfiles;posfiles];
mycellarray ={};

for i = 1:size(allfiles,1)
    disp(sprintf('Processing review %d out of %d', i, size(allfiles,1)));
    myfile = allfiles{i};
    fid = fopen( myfile);
    s = textscan(fid,'%s','Delimiter','\n');
    mystr = '';
    for mycellindex = 1:size(s{1,1},1)
        mystr = strcat(mystr, s{1,1}{mycellindex});
    end
    fclose(fid);
    mycellarray{end+1} = mystr;
end
mycellarray = mycellarray';
%mycellarray2 = mycellarray;
%mycellarray = mycellarray(1:5);

n=70
tic;

featureVector = featurize(mycellarray, n, 0, 0);
featureVectorOrig = featureVector;
save('featureVectorn70.dump','featureVector')
toc;
disp('Naive Bayes - Multinomial');

%10 fold random permutation
Fresults = [];

for i = 1:10
    randomindices = randperm(2000);
    randomindices = randomindices(1:1800);
    otherindices = (1:2000)';
    testsetindex = setdiff(otherindices,randomindices)';
    trainingsetindex = randomindices ;
    trainingset = featureVector(trainingsetindex,:);
    traininglabel = labels(trainingsetindex,:);

    testset = featureVector(testsetindex,:);
    testlabel = labels(testsetindex,:);
    O1 = NaiveBayes.fit(trainingset,traininglabel,'dist','mn'); % or  'mvmn'
    C2 = O1.predict(testset);
    cMat2 = confusionmat(testlabel,C2);
    Fresults = [Fresults,F1measureConfusionMatrix(cMat2)];
end
%disp(n);
fprintf('F1-measure for Naive Bayes classifier = %0.5f\n', mean(Fresults))

toc;

tic;

disp('Naive Bayes - Bernoulli');
%10 fold validation
Fresults = [];
featureVector = bernoulli(featureVectorOrig);

for i = 1:10
    randomindices = randperm(2000);
    randomindices = randomindices(1:1800);
    otherindices = (1:2000)';
    testsetindex = setdiff(otherindices,randomindices)';
    trainingsetindex = randomindices ;
    trainingset = featureVector(trainingsetindex,:);
    traininglabel = labels(trainingsetindex,:);

    testset = featureVector(testsetindex,:);
    testlabel = labels(testsetindex,:);
    O1 = NaiveBayes.fit(trainingset,traininglabel,'dist','mn'); % or  'mvmn'
    C2 = O1.predict(testset);
    cMat2 = confusionmat(testlabel,C2);
    Fresults = [Fresults,F1measureConfusionMatrix(cMat2)];
end
%disp(n);
fprintf('F1-measure for Naive Bayes classifier = %0.5f\n', mean(Fresults))

toc;

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function handle = imageVisualise(imageVals, imageSize, transpose, negative, ...
				 scale)

% IMAGEVISUALISE Helper code for showing an image during 2-D visualisation.
% FORMAT
% DESC is a helper function for plotting image data using latent
% variable models.
% ARG imageValues : the values to set the image data to.
% ARG imageSize : the size of the image.
% ARG transpose : whether the resized image needs to be transposed
% (default 1, which is yes).
% ARG negative : whether to display the negative of the image
% (default 0, which is no).
% ARG scale : whether or not to use the imagesc function (defaults
% to 1, which is yes).
% RETURN handle : a the handle to the image data.
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004, 2006
%
% SEEALSO : imageModify, lvmResultsDynamic

% MLTOOLS

if nargin < 3
  transpose = 1;
end
if nargin< 4
  negative = 0;
end
if nargin < 5
  scale = 1;
end
if negative
  imageVals = -imageVals;
end
imageData = reshape(imageVals, imageSize(1), imageSize(2));
if transpose
  imageData = imageData';
end
if scale
  handle = imagesc(imageData);
else
  handle = image(imageData);
end
colormap gray

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ Q ] = vec2Qr( vec )
phi = -vec(1:3);
rho = -vec(4:6);

ph = norm(phi);
ph2 = ph*ph;
ph3 = ph2*ph;
ph4 = ph3*ph;
ph5 = ph4*ph;

cph = cos(ph);
sph = sin(ph);

rx = hat(rho);
px = hat(phi);

t1 = 0.5 * rx;
t2 = ((ph - sph)/ph3) * (px*rx + rx*px + px*rx*px);
m3 = ( 1 - 0.5 * ph2 - cph ) / ph4;
t3 = -m3 * (px*px*rx + rx*px*px - 3*px*rx*px);
m4 = 0.5 * ( m3 - 3*(ph - sph - ph3/6)/ph5 );
t4 = -m4 * (px*rx*px*px + px*px*rx*px);

Q = t1 + t2 + t3 + t4;

end

author: mbrossar


text:
function c = sigma_delta_control

c = hdlnewcontrol(mfilename);

c.forEach('design_multi_stage3/Stage1_40_Tap_FIR_using_LUT/Sum',...
 'built-in/Sum', {},...
 'hdldefaults.ProductTreeHDLEmission',{});

%for tree implementation use hdldefaults.ProductTreeHDLEmission
%for cascade implementaiton use hdldefaults.ProductCascadeHDLEmission


text:
function c = matRad_constraintFunctions(optiProb,w,dij,cst)
% matRad IPOPT callback: constraint function for inverse planning
% supporting max dose constraint, min dose constraint, min mean dose constraint,
% max mean dose constraint, min EUD constraint, max EUD constraint,
% max DVH constraint, min DVH constraint
%
% call
%   c = matRad_constraintFunctions(optiProb,w,dij,cst)
%
% input
%   optiProb:   option struct defining the type of optimization
%   w:          bixel weight vector
%   dij:        dose influence matrix
%   cst:        matRad cst struct
%
% output
%   c:          value of constraints
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2016 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get current dose / effect / RBExDose vector
%d = matRad_backProjection(w,dij,optiProb);
optiProb.BP = optiProb.BP.compute(dij,w);
d = optiProb.BP.GetResult();

% Initializes constraints
c = [];

% compute objective function for every VOI.
for  i = 1:size(cst,1)

    % Only take OAR or target VOI.
    if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') )

        % loop over the number of constraints for the current VOI
        for j = 1:numel(cst{i,6})

            obj = cst{i,6}{j};

            % only perform computations for constraints
            % if ~isempty(strfind(obj.type,'constraint'))
            if isa(obj,'DoseConstraints.matRad_DoseConstraint')

                % if we have effect optimization, temporarily replace doses with effect
                % Maybe we should put some switch into the classes for that
                if (~isequal(obj.name, 'max dose constraint')      && ~isequal(obj.name, 'min dose constraint')      &&...
                    ~isequal(obj.name, 'max mean dose constraint') && ~isequal(obj.name, 'min mean dose constraint') && ...
                    ~isequal(obj.name, 'min EUD constraint')       && ~isequal(obj.name, 'max EUD constraint'))     && ...
                    (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection'))

                    doses = obj.getDoseParameters();

                    effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2;

                    obj = obj.setDoseParameters(effect);
                end

                % if conventional opt: just add constraints of nominal dose
                %if strcmp(cst{i,6}(j).robustness,'none')

                    d_i = d{1}(cst{i,4}{1});

                    %c = [c; matRad_constFunc(d_i,cst{i,6}(j),d_ref)];
                    c = [c; obj.computeDoseConstraintFunction(d_i)];

                %else

                    %error('Invalid robustness setting.');

                %end % if we are in the nominal sceario or rob opt

            end % if we are a constraint

        end % over all defined constraints & objectives

    end % if structure not empty and oar or target

end % over all structures

repo: matRad


text:
classdef IsotropicElasticMaterial < ElasticMaterial

    properties (GetAccess = public, SetAccess = protected)
        nstre
    end

    properties (Access = protected)
        kappa
        mu
        lambda
    end

    methods (Access = public)

        function compute(obj,s)
            obj.kappa = s.kappa;
            obj.mu    = s.mu;
            obj.nElem = size(obj.mu,1);
            obj.nGaus = size(obj.mu,2);
            obj.computeC();
        end

    end

    methods (Access = protected)

        function init(obj,cParams)
            obj.nstre = cParams.nstre;
        end

    end

    methods (Access = protected, Abstract)
        computeC(obj)
    end

    methods (Access = public, Static)

        function mu = computeMuFromYoungAndNu(E,nu)
            mu = E./(2*(1+nu));
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
 function [alphas, beta] = nufft_alpha_kb_fit(N, J, K, L, beta, chat)
%
% return the alpha and beta corresponding to LS fit of L components
% to optimized Kaiser-Bessel scaling factors (m=0, alpha=2.34J).
% This is the best method I know currently for choosing alpha!
%
% Copyright 2002-7-16, Jeff Fessler, The University of Michigan

if nargin < 3, help(mfilename), error(mfilename), end
if ~isvar('L') | isempty(L)
	if N > 40
		L = 13;		% empirically found to be reasonable
	else
		L = ceil(N/3);	% a kludge to avoid "rank deficient" complaints
	end
end
if ~isvar('beta') | isempty(beta),	beta = 1; end
if ~isvar('chat') | isempty(chat),	chat = 0; end

kb_alf = 2.34 * J;	% KB shape parameter
kb_m = 0;		% KB order

[tmp, sn_kaiser] = nufft1_error(0, N, J, K, 'kaiser', 'ft');
sn_kaiser = reale(sn_kaiser);

%
% use regression to match NUFFT with BEST kaiser scaling's
%
gam = 2*pi/K;
nlist = [0:(N-1)]' - (N-1)/2;
X = cos(beta*gam*nlist*[0:L]);	% [N,L]
coef = (X \ sn_kaiser)';	% regress(sn_kaiser, X)';
alphas = [reale(coef(1)) coef(2:end)/2];

if chat
	printf('cond # for LS fit to KB scale factors: %g', cond(X))
end

author: thomaskuestner


text:
% Studies the recovery probability of rank aware thresholding
% with sparsity level and number of channels in joint recovery.

close all;
clear all;
clc;
rng('default');
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');
% Signal space
N = 256;
% Number of measurements
M = 32;
% Number of signals
Ss = [1, 2, 4, 8, 16, 32];
% Sparsity levels
Ks = 2:30;
num_trials = 200;
num_ks = length(Ks);
num_ss = length(Ss);
success_with_k = zeros(num_ss, num_ks);

snr_threshold = 100;

for ns=1:num_ss
    % Current sparsity level
    S = Ss(ns);
    for nk=1:num_ks
        K = Ks(nk);
        num_successes = 0;
        for nt=1:num_trials
            % Sensing matrix
            Phi = spx.dict.simple.gaussian_dict(M, N);
            % Sparse signal generator
            gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
            % Gaussian distributed non-zero samples
            X = gen.gaussian;
            % Measurement vectors
            Y = Phi * X;
            % Create the solver for rank aware thresholding
            solver = spx.pursuit.joint.RankAwareThresholding(Phi, K);
            result = solver.solve(Y);
            % Solution vectors
            X_Rec = result.Z;
            % Comparison
            cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
            % Reconstruction SNR
            snr = cs.cum_signal_to_noise_ratio;
            success = snr > snr_threshold;
            num_successes = num_successes + success;
            fprintf('S: %d, K=%d, trial=%d, residual: %e, SNR: %f dB\n'...
                , S, K, nt, cs.cum_difference_norm, snr);
        end
        success_with_k(ns, nk) = num_successes / num_trials;
    end
end

save ('bin/figure_1_ra_thresholding_success_with_K_S.mat');

author: indigits

repo: sparse-plex


text:
function prtUtilClassStructurePlot(classNamesStrCell)
% prtUtilClassStructurePlot(classNamesStrCell);
%
% Mostly Ripped from viewClassTree() from MATLAB central by Matthew Dunham.
%
% View a class inheritence hierarchy. All classes residing in the directory
% or any subdirectory are discovered. Parents of these classes are also
% discovered as long as they are in the matlab search path.
% There are a few restrictions:
% (1) classes must be written using the new 2008a classdef syntax
% -This is no longer true-(2) classes must reside in their own @ directories.
% (3) requires the bioinformatics biograph class to display the tree.
% (4) works only on systems that support 'dir', i.e. windows.
%
%Written by Matthew Dunham

suppliedDir = false;
if nargin == 0
    directory = pwd;
    suppliedDir = true;
elseif ischar(classNamesStrCell)
    if isdir(classNamesStrCell)
        directory = classNamesStrCell;
        suppliedDir = true;
    else
        classNamesStrCell = {classNamesStrCell};
    end
else
    classNamesStrCell = {class(classNamesStrCell)};
    suppliedDir = false;
end

if suppliedDir
    %fileNames = what(directory);
    fileNames = prtUtilRecursiveWhat(directory);
    classNamesStrCell = cellfun(@(c)c(1:end-2), fileNames.m,'uniformOutput',false);
end

allClasses = classNamesStrCell;
for c=1:numel(classNamesStrCell)
    allClasses = union(allClasses,ancestors(classNamesStrCell{c}));
end
allClasses = unique(allClasses);

matrix = zeros(numel(allClasses));
map = struct;
for i=1:numel(allClasses)
    map.(allClasses{i}) = i;
end

keepClass = false(numel(allClasses),1);
for i=1:numel(allClasses)
    try
        m = eval(['?',allClasses{i}]);
        %m = meta(allClasses{i});
        parents = m.SuperClasses;
        if ~isempty(m)
            keepClass(i) = true;
        end
    catch ME %#ok
        warning('CLASSTREE:discoveryWarning',['Could not discover information about class ',allClasses{i}]);
        continue;
    end
    for j=1:numel(parents)
        matrix(map.(allClasses{i}),map.(parents{j}.Name)) = 1;
    end
end

matrix = matrix(keepClass,keepClass);
allClasses = allClasses(keepClass);

if isscalar(matrix)
    matrix = 1;
end

plotGraph(matrix, allClasses)

end

function info = dirinfo(directory)
%Recursively generate an array of structures holding information about each
%directory/subdirectory beginning, (and including) the initially specified
%parent directory.
info = what(directory);
flist = dir(directory);
dlist =  {flist([flist.isdir]).name};
for i=1:numel(dlist)
    dirname = dlist{i};
    if(~strcmp(dirname,'.') && ~strcmp(dirname,'..'))
        info = [info, dirinfo([directory,'\',dirname])]; %#ok
    end
end
end

function list = ancestors(className)
%Recursively generate a list of all of the superclasses, (and superclasses
%of superclasses, etc) of the specified class.
list = {};

try
    m = eval(['?',className]);
    parents = m.SuperclassList;
catch %#ok
    disp('error')
    return;
end

list = getAncestorsFromMeta(parents);

list = unique(list);

end

function list = getAncestorsFromMeta(m)
list = {};

for p = 1:numel(m)
    if(p > numel(m)),continue,end %bug fix for version 7.5.0 (2007b)
    list = cat(2, list, {m(p).Name}, getAncestorsFromMeta(m(p).SuperclassList));
end
end

function plotGraph(connectivity, nodeStrs)

% plotGraph a simple quick plot of the graph layout for validation.

% Call GraphViz to get a good layout - Call the MATLAB version
graphLayoutInfo = prtPlotUtilGraphVizRun(connectivity');

% Plot each node as a box with a string of it's actionId
nNodes = length(graphLayoutInfo.Nodes);
for iNode = 1:nNodes
    cN = graphLayoutInfo.Nodes(iNode);

    cNHandles.rect = rectangle('Position',[cN.x cN.y cN.w cN.h]);
    cNHandles.text = text(cN.x+cN.w/2,cN.y+cN.h/2, nodeStrs{iNode}, 'VerticalAlignment','Middle','HorizontalAlignment','Center','Interpreter','none');

    if iNode == 1
        nodeHandles = repmat(cNHandles,nNodes,1);
    else
        nodeHandles(iNode) = cNHandles;
    end
end

% Get the scale of the axes so we know how big to make arrow heads
v = axis;
scale = mean([v(2)-v(1),v(4)-v(3)]);

% Plot each edge
nEdges = length(graphLayoutInfo.Edges);
for iEdge = 1:nEdges
    cE = graphLayoutInfo.Edges(iEdge);

    cStartN = graphLayoutInfo.Nodes(cE.startIndex);
    cStopN = graphLayoutInfo.Nodes(cE.stopIndex);

    cArrowStart = cat(2,cStartN.x + cStartN.w/2, cStartN.y + cStartN.h);
    cArrowStop = cat(2,cStopN.x + cStopN.w/2, cStopN.y);

    cEHandles.arrow = prtPlotUtilPlotArrow(cat(1,cArrowStart(1),cArrowStop(1)), cat(1,cArrowStart(2),cArrowStop(2)), [],'headWidth',0.01*scale,'realHeadLength',0.01*scale);

    if iNode == 1
        edgeHandles = repmat(cEHandles,nEdges,1);
    else
        edgeHandles(iNode) = cEHandles;
    end
end

set(gca,'XTick',[],'YTick',[]);
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [s_train, s_test, vocab_item, vocab_user, ui_train] = load_data(train_path, test_path, L, T)
    % Codes for loading data from local file
    %
    % ARGS:
    % train_path    : the given path to training file
    % test_path     : the given path to test file
    % L             : the length of sequence
    % T             : number of targets
    %
    % RETURN:
    % s_train       : the resulted training instances
    % s_test        : the resulted test instances
    % vocab_item    : the relationship vocabulary between new item id and origin item id
    % vocab_user    : the relationship vocabulary between new user id and origin user id
    % ui_train      : rated item of each user

    max_length = L + T;
    train_ratio = L/max_length;

    fid_train = fopen(train_path);
    fid_test = fopen(test_path);

    % train
    data = textscan(fid_train, '%f\t%f\t%f');
    fclose(fid_train);
    train = cell2mat(data);
    clear data;

    vocab_user = unique(train(:, 1));
    vocab_item = unique(train(:, 2));

    s_train = cell(1e+6, 3); % data, label, user
    s_test = cell(1e+6, 3);
    ui_train = cell(length(vocab_user), 1);

    cnt = 1;
    for i=1:length(vocab_user)
        new_user_id = i;
        old_user_id = vocab_user(i);
        ind = train(:, 1) == old_user_id;
        items = train(ind, 2);
        [~, items_] = ismember(items, vocab_item);
        ui_train{new_user_id} = items_;
        if length(items_) <= max_length
            items_new = zeros(max_length,1);
            items_new(max_length - length(items_)+1:end) = items_;
            ind_train = floor(length(items_new) * train_ratio);

            s_train{cnt, 1} = items_new(1:ind_train);
            s_train{cnt, 2} = items_new(ind_train+1:end);
            s_train{cnt, 3} = new_user_id;
            cnt = cnt + 1;

            s_test{new_user_id, 1} = items_new(end-floor(max_length*train_ratio)+1:end);
            s_test{new_user_id, 3} = new_user_id;
        else
            ind_label = max_length:1:length(items_);
            for j=ind_label
                inds = (j-max_length+1) : j;
                ind_ = floor(max_length * train_ratio);
                ind_train = inds(1: ind_);
                ind_test = inds(ind_+1:end);
                s_train{cnt, 1} = items_(ind_train);
                s_train{cnt, 2} = items_(ind_test);
                s_train{cnt, 3} = new_user_id;
                cnt = cnt + 1;
            end

            s_test{new_user_id, 1} = items_(end-floor(max_length*train_ratio)+1:end);
            s_test{new_user_id, 3} = new_user_id;
        end
    end

    % test
    data = textscan(fid_test, '%f\t%f\t%f');
    fclose(fid_test);
    test = cell2mat(data);
    clear data;

    for i=1:length(vocab_user)
        new_user_id = i;
        old_user_id = vocab_user(i);

        ind = test(:, 1) == old_user_id;
        items = test(ind, 2);
        [~, items_] = ismember(items, vocab_item);
        s_test{new_user_id, 2} = items_;
    end

    empty_cells = cellfun(@isempty, s_train);
    empty_rows = empty_cells(:,1) == 1;
    s_train(empty_rows, :) = [];

    empty_cells = cellfun(@isempty, s_test);
    empty_rows = empty_cells(:,1) == 1;
    s_test(empty_rows, :) = [];
end

author: graytowne

repo: caser

save_path: github-repos/MATLAB/graytowne-caser


text:
function plot_sound_field(P,X,Y,Z,x0,conf)
%PLOT_SOUND_FIELD plots a sound field
%
%   Usage: plot_sound_field(P,X,Y,Z,[x0],conf)
%
%   Input parameters:
%       P           - matrix containing the sound field in the format P = P(y,x)
%       X           - x-axis / m; single value or [xmin,xmax] or nD-array
%       Y           - y-axis / m; single value or [ymin,ymax] or nD-array
%       Z           - z-axis / m; single value or [zmin,zmax] or nD-array
%       x,y,z       - vectors for the x-, y- and z-axis
%       x0          - matrix containing the secondary source positions to plot.
%                     Default: plot no secondary sources
%       conf        - configuration struct (see SFS_config)
%
%   PLOT_SOUND_FIELD(P,X,Y,Z,x0,conf) plots the sound field P in dependence of
%   the axes that are non-singleton. You have to provide the axes in the same
%   syntax as for the sound field functions. For a given set x0 of secondary
%   sources the secondary sources are added as dots or loudspeaker symbols
%   depending on your setting of conf.plot.realloudspeakers.
%
%   See also: sound_field_mono, sound_field_imp, draw_loudspeakers

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 6;
narginchk(nargmin,nargmax);
isargnumeric(X,Y,Z,P);
if nargin<nargmax
    conf = x0;
    x0 = [];
end
if isempty(x0)
    conf.plot.loudspeakers = 0;
end
isargstruct(conf);

%% ===== Configuration ==================================================
% Tmp dir
tmpdir = conf.tmpdir;
% Plotting
p.usenormalisation = conf.plot.usenormalisation;
p.cmd = conf.plot.cmd;
p.usedb = conf.plot.usedb;
p.mode = conf.plot.mode;
p.size = conf.plot.size;
p.size_unit = conf.plot.size_unit;
p.caxis = conf.plot.caxis;
p.colormap = conf.plot.colormap;
p.loudspeakers = conf.plot.loudspeakers;
p.realloudspeakers = conf.plot.realloudspeakers;
p.lssize = conf.plot.lssize;
p.usefile = conf.plot.usefile;
p.file = conf.plot.file;

%% ===== Calculation =====================================================
% Handle the given axis and check which should be plotted
[dimensions,x1,x2,x3] = xyz_axes_selection(X,Y,Z);
% Check whether some of the axis have custom size
% NOTE: the grid will always be detected as custom, if x,y,z have more than two
% entries, which could also happen if you perform for example:
% >> [P,x,y,z,x0] = sound_field_mono_wfs(X,Y,Z,[0 -1 0],'pw',1000,conf);
% >> plot_sound_field(P,x,y,z,x0,conf)
% Instead you should do
% >> plot_sound_field(P,X,Y,Z,x0,conf)
% if X,Y,Z are a non-custom grid.
is_custom = is_dim_custom(X,Y,Z);

if any(is_custom) && sum(is_custom) < sum(dimensions)
  error(['%s: it is not possible to combine an axis with regular sampled ' ...
    'grid with an axis with custom grid']);
end

if ~any(dimensions)
    error(['%s: you have only one point in the sound field. ', ...
        'Omitting the plotting.'],upper(mfilename));
elseif ~dimensions(1) && ~dimensions(2)
    p.xlabel = 'z / m';
elseif ~dimensions(1)
    p.xlabel = 'y / m';
    p.ylabel = 'z / m';
elseif ~dimensions(2)
    p.xlabel = 'x / m';
    p.ylabel = 'z / m';
elseif ~dimensions(3)
    p.xlabel = 'x / m';
    p.ylabel = 'y / m';
else
    p.xlabel = 'x / m';
    p.ylabel = 'y / m';
    p.zlabel = 'z / m';
end

% Remove Inf values from sound field
P(isinf(P)) = NaN;

% Normalisation
if p.usenormalisation
    P = norm_sound_field(P,conf);
end

if p.usedb
    % If we have only zeros in the sound field set the field to eps to avoid
    % problems with log(0).
    if max(abs(P(:)))==0
        P(:) = eps;
    end
    % Change default colormap to yellowred
    if strcmp('default',p.colormap)
        conf.plot.colormap = 'yellowred';
    end
    % Calculate sound pressure level in dB
    P = 20*log10(abs(P));
    if p.usenormalisation
        P = P - max(P(:)); % ensure 0 dB max for dB plots
    end
else
    P = real(P);
end

% Check if we should plot loudspeakers symbols
if p.realloudspeakers && size(x0,1)>1000
    warning(['%s: the given number of loudspeaker is >1000. ',...
            'Switching back to non loudspeaker symbols.'],upper(mfilename));
    p.realloudspeakers = 0;
end

% Set the color bar axis to default values if not given otherwise
if p.caxis, else
    if p.usedb
        p.caxis = [-45,0];
    else
        p.caxis = [-1,1];
    end
end

%% ===== Plotting ========================================================

% Create a new figure
figure;
% Set size
figsize(p.size(1),p.size(2),p.size_unit);

switch sum(dimensions)
  case 1
    % === 1D Plot ====
    if any(is_custom)  % custom grid
        plot(x1,P,'o','MarkerFaceColor','blue','MarkerSize',3)
    else % regular grid
        x1 = linspace(x1(1),x1(2),conf.resolution);
        plot(x1,P);
    end
    xlabel(p.xlabel);
    if p.usedb
        ylabel('Amplitude / dB');
    else
        ylabel('Amplitude');
    end
  case 2
    % === 2D Plot ====
    if any(is_custom)  % custom grid
        point_size = [];
        if isoctave, point_size = 2; end % fix scatter plot point size for Octave
        scatter(x1(:),x2(:),point_size,limit_colors(P,p.caxis),'filled');
    else  % regular grid
        imagesc(x1,x2,P,p.caxis);
    end
    % Add color bar
    set_colorbar(conf);
    % Set the y direction in normal mode (imagesc uses the reverse mode by
    % default)
    turn_imagesc;
    % Set the axis to use the same amount of space for the same length (m)
    axis image;
    % Labels etc. for the plot
    xlabel(p.xlabel);
    ylabel(p.ylabel);
    % Add loudspeaker to the plot
    if p.loudspeakers % && dimensions(1) && dimensions(2)
        hold on;
        draw_loudspeakers(x0,dimensions,conf);
        hold off;
    end
  case 3
    if ~any(is_custom)
      [x1,x2,x3] = xyz_grid(X,Y,Z,conf);
    end
    % === 3D Plot ====
    point_size = [];
    if isoctave, point_size = 5; end % fix scatter plot point size for Octave
    scatter3(x1(:),x2(:),x3(:),point_size,limit_colors(P,p.caxis),'filled');
    % Fix perspective of plot in Matlab (https://stackoverflow.com/q/17447404)
    set(gcf,'renderer','opengl');
    % Add color bar
    set_colorbar(conf);
    % Set the axis to use the same amount of space for the same length (m)
    axis image;
    % Labels etc. for the plot
    xlabel(p.xlabel);
    ylabel(p.ylabel);
    zlabel(p.zlabel);
end

% Save as file
if ~isempty(p.file) && strcmp('png',p.file(end-2:end))
    print_png(p.file);
elseif ~isempty(p.file) && strcmp('eps',p.file(end-2:end))
    print_eps(p.file);
end
end

%% ===== Helper functions ================================================
function P = limit_colors(P,caxis)
    % Limits the number of different values in P to 64.
    % This speeds up plotting with the scatter function in octave, see:
    % https://savannah.gnu.org/bugs/?40663
    %
    number_of_colors = 64;
    % First apply the caxis values
    P = min(caxis(2),max(caxis(1),P(:)));
    % Transform to [0...64] and transform to integer
    P = fix(number_of_colors/abs(max(P(:))-min(P(:))) .* (P+abs(min(P(:)))));
    % Transform back to [caxis(1)...caxis(2)]
    P = abs(caxis(2)-caxis(1))/number_of_colors .* P + caxis(1);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [no,el]=removeisolatednode(node,elem)
%
% [no,el]=removeisolatednode(node,elem)
%
% remove isolated nodes: nodes that are not included in any element
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
%
% input:
%     node: list of node coordinates
%     elem: list of elements of the mesh
%
% output:
%     no: node coordinates after removing the isolated nodes
%     el: element list of the resulting mesh
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

oid=1:size(node,1);       % old node index
idx=setdiff(oid,elem(:)); % indices to the isolated nodes
idx=sort(idx);
delta=zeros(size(oid));
delta(idx)=1;
delta=-cumsum(delta);     % calculate the new node index after removing the isolated nodes
oid=oid+delta;            % map to new index
el=oid(elem);             % element list in the new index
no=node;
no(idx,:)=[];             % remove the isolated nodes

repo: REST


text:
classdef TestFeatureDetector
    %TestFeatureDetector

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
        kfields = {'pt', 'size', 'angle', 'response', 'octave', 'class_id'};
        detectors = { ...
            ... % features2d (opencv)
            'AgastFeatureDetector', 'FastFeatureDetector', ...
            'GFTTDetector', 'SimpleBlobDetector', 'MSER', ...
            'BRISK', 'ORB', 'KAZE', 'AKAZE', ...
            ... % xfeatures2d (opencv_contrib)
            'SIFT', 'SURF', 'StarDetector', 'MSDDetector', ...
            'HarrisLaplaceFeatureDetector'
        };
    end

    methods (Static)
        function test_detect_img
            img = cv.imread(TestFeatureDetector.im, ...
                'Grayscale',true, 'ReduceScale',2);
            for i=1:numel(TestFeatureDetector.detectors)
                try
                    obj = cv.FeatureDetector(TestFeatureDetector.detectors{i});
                catch ME
                    %TODO: check if opencv_contrib/xfeatures2d is available
                    %error('mexopencv:testskip', 'contrib');
                    continue;
                end
                typename = obj.typeid();

                kpts = obj.detect(img);
                validateattributes(kpts, {'struct'}, {'vector'});
                assert(all(ismember(TestFeatureDetector.kfields, ...
                    fieldnames(kpts))));
            end
        end

        function test_detect_imgset
            img = cv.imread(TestFeatureDetector.im, ...
                'Grayscale',true, 'ReduceScale',2);
            imgs = {img, img};
            for i=1:numel(TestFeatureDetector.detectors)
                try
                    obj = cv.FeatureDetector(TestFeatureDetector.detectors{i});
                catch ME
                    %TODO: check if opencv_contrib/xfeatures2d is available
                    %error('mexopencv:testskip', 'contrib');
                    continue;
                end

                kpts = obj.detect(imgs);
                validateattributes(kpts, {'cell'}, ...
                    {'vector', 'numel',numel(imgs)});
                cellfun(@(kpt) validateattributes(kpt, ...
                    {'struct'}, {'vector'}), kpts);
                cellfun(@(kpt) assert(all(ismember(...
                    TestFeatureDetector.kfields, fieldnames(kpt)))), kpts);
            end
        end

        function test_detect_mask
            img = imread(TestFeatureDetector.im);
            mask = zeros(size(img), 'uint8');
            mask(:,1:end/2) = 255;  % only search left half of the image

            for i=1:numel(TestFeatureDetector.detectors)
                try
                    obj = cv.FeatureDetector(TestFeatureDetector.detectors{i});
                catch ME
                    %TODO: check if opencv_contrib/xfeatures2d is available
                    %error('mexopencv:testskip', 'contrib');
                    continue;
                end

                kpts = obj.detect(img, 'Mask',mask);
                validateattributes(kpts, {'struct'}, {'vector'});
                assert(all(ismember(TestFeatureDetector.kfields, ...
                    fieldnames(kpts))));
                xy = cat(1, kpts.pt);
                assert(all(xy(:,1) <= ceil(size(img,2)/2)));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function sparse_grid_mixed_unique_index_test ( dim_num, level_max_min, ...
  level_max_max, rule, alpha, beta, tol )

%*****************************************************************************80
%
%  Purpose:
%
%    SPARSE_GRID_MIXED_UNIQUE_INDEX_TEST tests SPARSE_GRID_MIXED_UNIQUE_INDEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX_MIN, LEVEL_MAX_MAX, the minimum and
%    maximum values of LEVEL_MAX.
%
%    Input, integer RULE(DIM_NUM), the rule in each dimension.
%     1, "CC",  Clenshaw Curtis, Closed Fully Nested rule.
%     2, "F2",  Fejer Type 2, Open Fully Nested rule.
%     3, "GP",  Gauss Patterson, Open Fully Nested rule.
%     4, "GL",  Gauss Legendre, Open Weakly Nested rule.
%     5, "GH",  Gauss Hermite, Open Weakly Nested rule.
%     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested rule.
%     7, "LG",  Gauss Laguerre, Open Non Nested rule.
%     8, "GLG", Generalized Gauss Laguerre, Open Non Nested rule.
%     9, "GJ",  Gauss Jacobi, Open Non Nested rule.
%    10, "GW",  Golub Welsch, (presumed) Open Non Nested rule.
%    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested rule.
%    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested rule.
%    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested rule.
%    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested rule.
%    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested rule.
%    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested rule.
%    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
%
%    Input, real ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
%    Generalized Gauss Hermite, Generalized Gauss Laguerre, and Gauss Jacobi rules.
%
%    Input, real TOL, the tolerance for point equality.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_UNIQUE_INDEX_TEST\n' );
  fprintf ( 1, '  SPARSE_GRID_MIXED_UNIQUE_INDEX returns a mapping between\n' );
  fprintf ( 1, '  the nonunique and unique points in a sparse grid.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, ' Dimension      Rule     Alpha          Beta\n' );
  fprintf ( 1, '\n' );

  for dim = 1 : dim_num
    fprintf ( 1, '  %8d  %8d', dim, rule(dim) );
    if ( rule(dim) == 6 | rule(dim) == 8 | rule(dim) == 9 )
      fprintf ( 1, '  %12e', alpha(dim) );
    end
    if ( rule(dim) == 9 )
      fprintf ( 1, '  %12e', beta(dim) );
    end
    fprintf ( 1, '\n' );
  end

  for level_max = level_max_min : level_max_max

    point_total_num = sparse_grid_mixed_size_total ( dim_num, level_max, rule );

    point_num = sparse_grid_mixed_size ( dim_num, level_max, rule, alpha, ...
      beta, tol );

    fprintf ( 1, '\n' );
    fprintf ( 1, ' LEVEL_MIN LEVEL_MAX POINT_NUM POINT_NUM\n' );
    fprintf ( 1, '                        Unique     Total\n' );
    fprintf ( 1, '\n' );

    level_min = max ( 0, level_max + 1 - dim_num );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %8d  %8d  %8d  %8d\n', ...
      level_min, level_max, point_num, point_total_num );

    sparse_unique_index = sparse_grid_mixed_unique_index ( ...
      dim_num, level_max, rule, alpha, beta, tol, point_num, point_total_num );

    fprintf ( 1, '\n' );
    fprintf ( 1, '     POINT    UNIQUE\n' );
    fprintf ( 1, '\n' );
    for point = 1 : point_total_num
      fprintf ( 1, '  %8d  %8d\n', point, sparse_unique_index(point) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function wavelenvsth
% _____________________________________________________
% This function gets the nomograph which shows  the relationship which exists
% between the wavelength, wave period and water depth using the dispersion
% equation.
%
% 0. Syntax:
% >> wavelenvsth
%
% 1. Inputs:
%     None.
%
% 2. Outputs:
%     Nomograph.
%
% 3. Example:
%  >> run(wavelenvsth)
%
%  4. Notes:
%     - It's only necessary to run the function.
%     - In the nomograph, it is activated the data cursor mode in order to
%       get the wavelength and waver depth in function of the water period.
%     - Wave depths in meters.
%     - Wave length in meters.
%
% 5. Referents:
%      Darlymple, R.G. and Dean R.A. (1999). Water Wave Mechanics
%              for Engineers and Scientist. Advanced Series on Ocean Engineering, Vol. 2
%              World Scientific. Singapure.
%      Le Mehuate, Bernard. (1976). An introduction to hidrodynamics and water waves.
%              Springer-Verlag. USA.
%
% Gabriel Ruiz.
% Jun-2006
% UNAM
%_______________________________________________________________________

%%%%%%%%%%  M A I N      F U N C T I O N   %%%%%%%%%%
% clear all; clc;
 h = 1:1:130;
    g = 9.81;
        T = 4.5;
            dh = h';
        m=1;
    fg =figure('Menubar', 'none', 'Name', 'Wavelength vs Period and Depth', 'NumberTitle', 'off',...
                    'Position' , [ 6 37 1011 697 ] , 'Color' , [ 0.87 0.87 0.87 ]);

    while T ~= 17

        for i=1:130
                    con = 1;
                        l(con) = 0;
                    l(con+1) = 1.56 * T ^ 2;

                    while abs( l(con+1) - l(con) ) > 0.0001,
                                l(con+2) = ( ( 9.81 * T ^ 2 ) / ( 2 * pi ) ) * tanh( ( 2 * pi * h(i) ) / l(con+1) );
                                con = con + 1;
                    end

                    L(i) =  l(con);
                    k = ( 2 * pi ) / L(i);
        end

        Ls{1,m} = L';
            hold on
                fgh =plot(dh,Ls{1,m});
                    xlabel('Water Depth (h), meters');
                        ylabel('Wave Length (L), meters');
                    xlim( [ 1 130 ] )
                set(fgh,'Color',rand(1,3));
            set(gca,'Box','on');
        drt = title('Wavelength vs Period and Depth', 'HorizontalAlignment' , 'center' , 'FontWeight', 'bold');
            set(gca, 'XGrid', 'off', 'XMinorTick', 'on' , 'YGrid' , 'off' , 'YMinorTick' , 'on', 'Fontsize', 8 );
                m = m+1;
                    T = 0.5+T;
    end

    T= 4.5:0.5:17;
        hj =length(T);
    textos2 = 'T = ';

    for i = 1 : hj
            stringer{i,1} = num2str(T(1,i));
            textos{i,1} = horzcat(textos2,stringer{i,1},' s');
    end
        asdk = legend(textos, 'Location', 'EastOutside');
            set(asdk, 'FontSize', 6);
                clc;
            datacursormode on;
        dcm_obj = datacursormode(fg);
       fundat = str2func('camdatos');
    set(dcm_obj, 'DisplayStyle', 'Window' , 'UpdateFcn' , fundat)
%%%%%%%%%%%%    E N D     M A I N     F U N C T I O N %%%%%%%%%%%%%

function [txtdcm,pos] = camdatos(empt,event_obj)    %% Subfunction 1
pos = get(event_obj,'Position');
txtdcm = {['Water Depth: ',num2str(pos(1))] , ['Wave Length: ',num2str(pos(2))]};

% End of Subfunction 1

author: Sable

repo: mcbench-benchmarks


text:
function [gx,dgdx] = g_odds(x,P,u,in)

x = x-max(x);
gx = (exp(x)+0)./sum(exp(x)+0);

n = size(x,1);
I = eye(n);
dgdx = (I - gx*ones(1,n))*diag(gx);

author: MBB-team

repo: VBA-toolbox


text:
3     % problem
3     % cavity type
7     % grid parameter
1     % uniform/stretched grid
1     % discretisation
0.2   % viscosity parameter
1     % Picard/Newton/hybrid linearization
9     % number of Picard iterations
1.d-5 % nonlinear tolerance
2     % uniform/exponential streamlines

%% Data file for driven cavity


text:
%                               431-400 Year Long Project
%                               LA1 - Medical Image Processing 2003
% Supervisor     :  Dr Lachlan Andrew
% Group Members  :  Alister Fong    78629   a.fong1@ugrad.unimelb.edu.au
%                   Lee Siew Teng   102519  s.lee1@ugrad.unimelb.edu.au
%                   Loh Jien Mei    103650  j.loh1@ugrad.unimelb.edu.au
%
% File and function name : calculate_min_distance
% Version                : 1.0
% Date of completion     : 6 October 2003
% Written by    :   Alister Fong    78629   a.fong1@ugrad.unimelb.edu.au
%
% Input   :
%           edge1,edge2 -   [X,Y] coordinates to be compared and measured
%           'testing' or 'not testing' (optional) - Default 'not testing'
%
% Output  :
%           min_distance - Returns the minimum distance between the two
%                          sets of coordinates.
%           matching_coordinates - [edge1X,edge1Y,edge2X,edge2Y]
%                           Pairs of coordinates of edge1 and edge 2 that
%                           are of the minimum distance.
%
% Description:
%       Calculates the minimum distance between 2 edges and returns the distance and
%   the matrix linking the coordinates to both edge coordinates.
%
% Usage >> [min_distance,matching_coordinates] = calculate_min_distance(edge1,edge2)
%                   or
%          [min_distance,matching_coordinates] = calculate_min_distance(edge1,edge2,'testing')
%                   or
%          [min_distance,matching_coordinates] = calculate_min_distance(edge1,edge2,'not testing')
%
% Example >> [min_distance,matching_coordinates] = calculate_min_distance(edge1,edge2)
%                 figure;
%                 plot(edge1(:,1),edge1(:,2),'r+-');
%                 hold on;
%                 plot(edge2(:,1),edge2(:,2),'b+-');
%                 for n = 1:1:length(matching_coordinates(:,1))
%                     plot([matching_coordinates(n,1);matching_coordinates(n,3)],...
%                          [matching_coordinates(n,2);matching_coordinates(n,4)],'g*-');
%                 end
%                 xlabel(strcat('minimum distance = ',num2str(min_distance)));
%
% WARNING >> In some coordinates there may visually be some overlaping lines that do not register
%            as being of distance 0. This is because the pixels representing these lines
%            do not coincide.

function [min_distance,matching_coordinates] = calculate_min_distance(edge1,edge2,varargin)
% ---------------------------------------------------------------------------------
% Process the input
% ---------------------------------------------------------------------------------
testing = 'not testing';
if ~isempty(varargin)
    for n = 1:1:length(varargin)
        if strcmp(varargin{n},'not testing') | strcmp(varargin{n},'testing')
            testing = varargin{n};
        end
    end
end
if isempty(edge1) | isempty(edge2)
    error('Input entry is a null matrix');
end
% ---------------------------------------------------------------------------------
% Find the closest edge and calculate the distance
% ---------------------------------------------------------------------------------
temp1 = [];
temp2 = [];
for n = 1:1:length(edge1(:,1))
    temp = edge2;
    temp(:,1) = edge1(n,1);
    temp(:,2) = edge1(n,2);
	temp1 = [temp1;temp];
    temp2 = [temp2;edge2];
    % Distance calculated here to reduce memory usage
    distance = euclidean_distance(temp1,temp2);
    min_distance = min(distance);
    pos = find(distance == min_distance);
    temp1 = temp1(pos,:);
    temp2 = temp2(pos,:);
end
distance = euclidean_distance(temp1,temp2);
pos = find(distance == min_distance);
min_distance = min(distance);

matching_coordinates = [temp1(pos,:),temp2(pos,:)];

% ---------------------------------------------------------------------------------
% Displaying the results for testing
% ---------------------------------------------------------------------------------
if strcmp(testing,'testing')
    figure;
    plot(edge1(:,1),edge1(:,2),'r+');
    hold on;
    plot(edge2(:,1),edge2(:,2),'b+');
    for n = 1:1:length(matching_coordinates(:,1))
        plot([matching_coordinates(n,1);matching_coordinates(n,3)],...
             [matching_coordinates(n,2);matching_coordinates(n,4)],'g*-');
    end
    title('edge1(red+), edge2(blue+) and selected edges(green*)');
    xlabel(strcat('minimum distance = ',num2str(min_distance)));
end


text:
function y = logsumexp( varargin )

%LOGSUMEXP    log(sum(exp(x))).
%   LOGSUMEXP(X) = LOG_SUM_EXP(X) = LOG(SUM(EXP(X)). We have replaced this
%   function with LOG_SUM_EXP to better match our function naming
%   conventions. Please start using it instead.

warning( 'CVX:Renamed', [ ...
    'The function "logsumexp" has been renamed "log_sum_exp". Please start\n', ...
    'using the new name. The old name will be removed in a future release.' ], 1 );

y = log_sum_exp( varargin{:} );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
function d = bpdhe(im)

im = im2uint8(im);

hsv = rgb2hsv(im);
h = hsv(:,:,1);
s = hsv(:,:,2);
i = uint8(hsv(:,:,3).*255);
% i = uint8(20+rand(512).*(200));
ma = double(max(i(:)));
mi = double(min(i(:)));
bins = (ma-mi)+1;
hist_i = hist(double(i(:)),bins);
gausFilter = fspecial('gaussian',[1 9],1.0762);
blur_hist = (imfilter(hist_i,gausFilter,'replicate'));
derivFilter = [-1 1];
deriv_hist = imfilter(blur_hist,derivFilter,'replicate');
sign_hist = sign(deriv_hist);
meanFilter = [1/3 1/3 1/3];
smooth_sign_hist = sign(imfilter(sign_hist,meanFilter,'replicate'));
cmpFilter = [1 1 1  -1 -1 -1 -1 -1];
index = zeros([1,3]);
index(1) = 0;
p = 2;
for n = 1:bins-7
    C = smooth_sign_hist(n:n+7) == cmpFilter;
    if sum(C) ==8
        index(p) = n+3;
        p = p+1;
    end
end
index(p) = bins;
factor = zeros([length(index)-1,1]);
span = factor;
M = factor;
range = factor;
start = factor;
endd = factor;
sub_hist = cell([length(index)-1,1]);
for m = 1:length(index)-1;
    sub_hist{m} = hist_i(index(m)+1:index(m+1));
    M(m) = sum(sub_hist{m});
    low = mi + index(m);
    high = mi + index(m+1)-1;
    span(m) = high-low+1;
    factor(m) = span(m)*log10(M(m));
end
factor_sum = sum(factor);
for m = 1:length(index)-1;
    range(m) = round((256-mi)*factor(m)/factor_sum);
end
start(1) = mi;
endd(1) = mi+range(1)-1;
for m = 2:length(index)-1;
    start(m) = start(m-1)+range(m-1);
    endd(m) = endd(m-1)+range(m);
end
y = cell([length(index)-1,1]);
s_r = zeros([1,mi]);
for m = 1:length(index)-1;
    hist_cum = cumsum(sub_hist{m});
    c = hist_cum./M(m);
    y{m} = round(start(m)+(endd(m)-start(m)).*c);
    s_r = [s_r,y{m}];
end
i_s = zeros(size(i));
for n = mi:ma
    lc = i == n;
    i_s(lc) = double(s_r(n+1))/255;
end
% hist_is = hist(double(i_s(:)),bins);
hsi_o = cat(3,h,s,i_s);
d = uint8(hsv2rgb(hsi_o).*255);

author: dawnlh

repo: awesome-low-light-image-enhancement


text:
function P_density = getInputPowerDensity(t,t0,tf,x,param,extra)
%	getInputPowerDensity returns the value of the input current density as a function of time.
%
%       P = getInputPowerDensity(t,t0,tf,extra)
%
%       Inputs:
%               - t     : value of the current time step
%               - t0    : initial integration time
%               - tf    : final integration time
%               - extra : extra parameters
%       Outputs:
%               - P_density     : Applied power density [W/m^2]

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Define your own linear/nonlinear function of time for the applied power density

% P = (t-t0)/(tf-t0) *(-30) + 0;
% P_density = -90*sin(t+100/100);

% P_density = -10*sin(omega*t+phi);
frequency = 0.1;
P_density = -300*sin((2*pi*frequency)*t+(1/100));
if P_density == 0
    P_density = 1e-2;
end

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function y = rz_func(x)

% nodal crossing function

% input

%  x = simulation time argument (days)

% output

%  y = z-component of unit position vector

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global jdtdbi iplanet1

% compute planet's state vector

jdtdb = jdtdbi + x;

[r, v] = pecliptic(jdtdb, iplanet1, 11);

% z-component of unit position vector

y = r(3) / norm(r);