text:
function g = linKernGradient(kern, x, varargin)

% LINKERNGRADIENT Gradient of LIN kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% linear
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO linKernParamInit, kernGradient, linKernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2009

% KERN

if nargin < 4
  [k, sk] = linKernCompute(kern, x);
else
  [k, sk] = linKernCompute(kern, x, varargin{1});
end
g(1) = sum(sum(varargin{end}.*sk));


text:
function errMsg = checkIfInputBodiesAreSame(departBody, arriveBody, prevErrMsgs, varargin)
%checkIfInputBodiesAreSame Summary of this function goes here
%   Detailed explanation goes here
    if(nargin==5)
        bodyName1 = varargin{1};
        bodyName2 = varargin{2};
    else
        bodyName1 = 'Departure body';
        bodyName2 = 'arrival body';
    end
    errMsg = prevErrMsgs;
    if(strcmpi(departBody, arriveBody))
        departBody = ['(', cap1stLetter(lower(departBody)), ')'];
        arriveBody = ['(', cap1stLetter(lower(arriveBody)), ')'];
        errMsg{end+1} = [bodyName1, ' ', departBody, ' and ', bodyName2, ' ', arriveBody, ' must be different.'];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function definput=arg_thresh(definput)

  definput.flags.iofun={'hard','soft','wiener'};
  definput.flags.outclass={'full','sparse'};

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [px,py,pz,tx,ty,tz,bx,by,bz,ierr] = nd2pt(wanx,wany,wanz,wdx,wdy,wdz)

%% c     compute Cartesian component of P, T and B axes from outward normal
%% c     and slip vectors
%% c
%% c     usage:
%% c     call nd2pt(anx,any,anz,dx,dy,dz,px,py,pz,tx,ty,tz,bx,by,bz,ierr)
%% c
%% c     arguments:
%% c     anx,any,anz    components of fault plane outward normal vector in the
%% c                    Aki-Richards Cartesian coordinate system (INPUT)
%% c     dx,dy,dz       components of slip vector in the Aki-Richards
%% c                    Cartesian coordinate system (INPUT)
%% c     px,py,pz       components of downward P (maximum dilatation) axis versor
%% c                    in the Aki-Richards Cartesian coordinate system (OUTPUT)
%% c     tx,ty,tz       components of downward T (maximum tension) axis versor
%% c                    in the Aki-Richards Cartesian coordinate system (OUTPUT)
%% c     bx,by,bz       components of downward B (neutral) axis versor in the
%% c                    Aki-Richards Cartesian coordinate system (OUTPUT)
%% c     ierr           error indicator (OUTPUT)
%% c
%% c     errors:
%% c     1              input vectors not perpendicular among each other

%% c
%% c      implicit none
%% c-------------------------------------------------------------------------------
%%       integer io
%%       real amistr,amastr,amidip,amadip,amirak,amarak,amitre,amatre
%%      1,amiplu,amaplu,orttol,ovrtol,tentol,dtor,c360,c90,c0,c1,c2,c3
%%       common /fpscom/amistr,amastr,amidip,amadip,amirak,amarak,amitre
%%      1,amatre,amiplu,amaplu,orttol,ovrtol,tentol,dtor,c360,c90,c0,c1,c2
%%      2,c3,io
%% c-------------------------------------------------------------------------------
%%       real wanx,wany,wanz,amn,anx,any,anz,wdx,wdy,wdz,amd,dx,dy,dz
%%      1,ang,px,py,pz,tx,ty,tz,bx,by,bz,amp
%%       integer ierr
%% c

%%       call fpsset
         amistr=-360.;
         amastr=360.;
         amidip=0.;
         amadip=90.;
         amirak=-360.;
         amarak=360.;
         amitre=-360.;
         amatre=360.;
         amiplu=0.;
         amaplu=90.;
         orttol=2.;
         ovrtol=0.001;
         tentol=0.0001;
         dtor=0.017453292519943296;
         c360=360.;
         c90=90.;
         c0=0.;
         c1=1.;
         c2=2.;
         c3=3.;
         io=6;

      ierr=0;
      [amn,anx,any,anz] = focal_norm(wanx,wany,wanz);
      [amd,dx,dy,dz] = focal_norm(wdx,wdy,wdz);
      [ang] = focal_angle(anx,any,anz,dx,dy,dz);
      if (abs(ang-c90) > orttol)
        disp(['ND2PT: input vectors not perpendicular, angle=' num2str(ang)]);
        ierr=1;
      end
      px=anx-dx;
      py=any-dy;
      pz=anz-dz;
      [amp,px,py,pz] = focal_norm(px,py,pz);
      if (pz < c0)
        [px,py,pz] = focal_invert(px,py,pz);
      end
      tx=anx+dx;
      ty=any+dy;
      tz=anz+dz;
      [amp,tx,ty,tz] = focal_norm(tx,ty,tz);
      if (tz < c0)
        [tx,ty,tz] = focal_invert(tx,ty,tz);
      end
      [bx,by,bz] = focal_vecpro(px,py,pz,tx,ty,tz);
      if(bz < c0)
        [bx,by,bz] = focal_invert(bx,by,bz);
      end

author: CelsoReyes


text:
function r = isempty(a)
%ISEMPTY      Returns 1 if input is empty, i.e. [], in the sense of Matlab
%
%   r = isempty(a)
%
%There are no empty intervals in the mathematical sense in INTLAB. For details,
%  please see Readme.txt
%

% written  09/29/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 10/12/05     S.M. Rump  functionality exchanged with @long\isempty_
% modified 04/22/09     S.M. Rump  isempty_ removed
%

  r = isempty(a.sign) | isempty(a.mantissa) | isempty(a.exponent);

repo: OpenMAS


text:
function approx_display ( )

%*****************************************************************************80
%
%% APPROX_DISPLAY displays a sequence of Bernstein approximants to SIN(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 July 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'APPROX_DISPLAY\n' );
  fprintf ( 1, '  BPAB_APPROX evaluates the Bernstein polynomial\n' );
  fprintf ( 1, '  approximant to a function F(X).\n' );
  fprintf ( 1, '  This program displays the sequence of approximants\n' );
  fprintf ( 1, '  to y=sin(x) over the interval [1,3]\n' );

  a = 1.0;
  b = 3.0;

  for degree = 0 : 20
%
%  Generate data values.
%
    xdata = zeros ( degree + 1, 1 );
    ydata = zeros ( degree + 1, 1 );

    for i = 0 : degree

      if ( degree == 0 )
        xdata(i+1) = 0.5 * ( a + b );
      else
        xdata(i+1) = ( ( degree - i ) * a   ...
                     + (          i ) * b ) ...
                     / ( degree     );
      end

      ydata(i+1) = sin ( xdata(i+1) );

    end
%
%  Compare the true function and the approximant.
%
    nval = 501;

    xval = linspace ( a, b, nval );

    yval = bpab_approx ( degree, a, b, ydata, nval, xval );

    clf
    plot ( xval, yval, 'b-' );
    hold on
    plot ( xval, sin ( xval ), 'r-' );
    plot ( xdata, ydata, 'r.', 'Markersize', 30 );
    grid
    xlabel ( '<---X--->' );
    ylabel ( '<---Y--->' );
    title ( sprintf ( 'Bernstein Approximant Degree %d', degree ) );
    hold off

    fprintf ( 1, 'Bernstein approximant of degree %d\n', degree );

    pause

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function S = radix4FFT1_Float(s)
    % This is a radix-4 FFT, using decimation in frequency
    % The input signal must be floating point
    % Works with real or complex input

    % Initialize variables and signals
    % NOTE: The length of the input signal should be a power of 4: 4, 16, 64, 256, etc.
    N = length(s);
    M = log2(N)/2;

    % Initialize variables for floating point sim
    W=exp(-j*2*pi*(0:N-1)/N);
    S = complex(zeros(1,N));
    sTemp = complex(zeros(1,N));

    % FFT algorithm
    % Calculate butterflies for first M-1 stages
    sTemp = s;
    for stage = 0:M-2
        for n=1:N/4
            S((1:4)+(n-1)*4) = radix4bfly(sTemp(n:N/4:end), floor((n-1)/(4^stage)) *(4^stage), 1, W);
        end
        sTemp = S;
    end

    % Calculate butterflies for last stage
    for n=1:N/4
        S((1:4)+(n-1)*4) = radix4bfly(sTemp(n:N/4:end), floor((n-1)/(4^stage)) * (4^stage), 0, W);
    end
    sTemp = S;

    % Rescale the final output
    S = S*N;

end

function Z = radix4bfly(x,segment,stageFlag,W)
    % For the last stage of a radix-4 FFT all the ABCD multiplers are 1.
    % Use the stageFlag variable to indicate the last stage
    % stageFlag = 0 indicates last FFT stage, set to 1 otherwise

    % Initialize variables and scale to 1/4
    a=x(1)*.25;b=x(2)*.25;c=x(3)*.25;d=x(4)*.25;

    % Radix-4 Algorithm
    A=a+b+c+d;
    B=(a-b+c-d)*W(2*segment*stageFlag + 1);
    C=(a-b*j-c+d*j)*W(segment*stageFlag + 1);
    D=(a+b*j-c-d*j)*W(3*segment*stageFlag + 1);

    Z = [A B C D];

end

author: Sable

repo: mcbench-benchmarks


text:
function [f] = elec_sphere_fit_optim(r, X, Y, Z, xo, yo, zo)

% elec_sphere_fit_optim - Optimization for elec_sphere_fit.m
%
% Called from elec_sphere_fit.m
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence:  GNU GPL, no implied or express warranties
% History:  02/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% with center (Xo,Yo,Zo) and radius r, the equation of a sphere is:
%
% r^2 = (x-xo)^2  +  (y-yo)^2  +  (z-zo)^2
%
% This function below creates a scalar value to
% return to the fminsearch function in elec_sphere_fit.

S = (X-xo).^2  +  (Y-yo).^2  +  (Z-zo).^2  -  r^2;

f = sum( S.^2 );

author: PatternRecognition

repo: OpenBMI


text:
function [ f ] = pendulumF( y )
%pendulumf Pendulum equations.

    g = 1;
    m = 1;
    l = 1;

    f(1,1) = y(3,:);
    f(2,1) = y(4,:);
    f(3,1) = -y(1,:).*y(5,:)/m;
    f(4,1) = (-y(2,:).*y(5,:)-g)/m;
    f(5,1) = y(1,:).*y(3,:) + y(2,:).*y(4,:);

end


text:
% convnet_classify
% Copyright (C) 2013 KyungHyun Cho
%
%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [c, posterior] = convnet_classify(C, x0, raw)

if nargin < 3
    raw = 0;
end

layers = C.structure.layers;
n_layers = length(layers);

posterior = x0;

full_layers = C.structure.full_layers;
conv_layers = C.structure.conv_layers;
n_full = length(full_layers);
n_conv = size(conv_layers,1);

strides = C.structure.strides;
poolratios = C.structure.poolratios;

mb_sz = size(x0, 1);

cin = C.structure.channel_in;
szin = C.structure.size_in;

repost = reshape(posterior, [mb_sz, szin, szin, cin]);

if C.lcn.use
    subwindow = fspecial('gaussian', C.lcn.neigh);
    %subwindow_sum = ones(C.lcn.neigh);
end

if C.lcn.use
    subsum = convn(repost, reshape(subwindow_sum, [1, C.lcn.neigh, C.lcn.neigh, 1]), 'same');
    repost = repost - subsum / C.lcn.neigh^2;
    repost2 = repost.^2;
    subsum = convn(repost2, reshape(subwindow_sum, [1, C.lcn.neigh, C.lcn.neigh, 1]), 'same');
    repost = repost ./(sqrt(subsum + 1e-12) / C.lcn.neigh);
end

for l = 1:n_conv
    fsz = length(C.cW{l}(:,1)) / cin;
    n_filters = size(C.cW{l}, 2);
    szout_prepool = (floor((szin - sqrt(fsz)) / strides(l)) + 1);
    szout = ceil(szout_prepool / poolratios(l));
    respfull = zeros(mb_sz, szout, szout, n_filters);

    % for each filter map
    for fidx = 1:n_filters
        fil = reshape(C.cW{l}(:, fidx), [1, sqrt(fsz), sqrt(fsz), cin]);
        resp = convn(repost, fil, 'valid');

        % stride
        resp = resp(:, 1:strides(l):end, :);
        resp = resp(:, :, 1:strides(l):end);

        resp = resp + C.cbiases{l}(fidx);

        % nonlinearity
        resp = sigmoid(resp, C.hidden.use_tanh);

        if C.lcn.use
            subsum = convn(resp, reshape(subwindow_sum, [1, C.lcn.neigh, C.lcn.neigh, 1]), 'same');
            resp = resp - subsum / C.lcn.neigh^2;
            resp2 = resp.^2;
            subsum = convn(resp2, reshape(subwindow_sum, [1, C.lcn.neigh, C.lcn.neigh, 1]), 'same');
            resp = resp ./(sqrt(subsum + 1e-12) / C.lcn.neigh);
        end

        if C.structure.poolratios(l) > 1
            % pooling
            switch C.pooling(l)
                case 0
                    resp = convnet_maxpool (resp, C.structure.poolratios(l));
                case 1
                    resp = convnet_avgpool (resp, C.structure.poolratios(l));
                case 2
                    error('NOT SUPPORTED');
            end
        end

        respfull(:,:,:,fidx) = resp;
    end

    szin = szout;
    cin = n_filters;
    repost = respfull;
end

% stretch
posterior = reshape(repost, [mb_sz C.structure.layers(n_conv+1)]);

for l = 2:n_full+1
    if C.dropout.use && l > 2
        posterior = (1 - C.dropout.prob) * posterior;
    end

    posterior = bsxfun(@plus, posterior * C.W{l-1}, C.biases{l}');

    if l < n_full + 1
        posterior = sigmoid(posterior, C.hidden.use_tanh);
    end
    if l == n_full + 1 && C.output.binary
        posterior = softmax(posterior);
    end
end

if raw
    c = posterior;
else
    [maxp, c] = max(posterior, [], 2);
end

author: kyunghyuncho

repo: deepmat


text:
function [vtx,elt,col] = mshReadMesh(filename)
%+========================================================================+
%|                                                                        |
%|                 OPENMSH - LIBRARY FOR MESH MANAGEMENT                  |
%|           openMsh is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : mshReadMesh.m                                 |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Read .mesh files (vertex, triangular elements,|
%|  `---'  |                and colours)                                  |
%+========================================================================+

% Open file
fid = fopen(filename,'r');
if( fid==-1 )
    error('Can''t open the file.');
end

% Nodes
str = fgets(fid);
while isempty(strfind(str,'Vertices'))
    str = fgets(fid);
end
Nvtx = str2double(fgets(fid));
vtx  = zeros(Nvtx,3);
for i = 1:Nvtx
    tmp = str2num(fgets(fid));
    vtx(i,:) = tmp(1:3);
end

% Support only triangular elements
str = fgets(fid);
if strcmp(str(1:9),'Triangles')
    Nelt = str2double(fgets(fid));
    elt  = zeros(Nelt,4);
    for i = 1:Nelt
        elt(i,:) = str2num(fgets(fid));
    end
    col = elt(:,4);
    elt = elt(:,1:3);

elseif strcmp(str(1:10),'Tetrahedra')
    Nelt = str2double(fgets(fid));
    elt  = zeros(Nelt,5);
    for i = 1:Nelt
        elt(i,:) = str2num(fgets(fid));
    end
    col = elt(:,5);
    elt = elt(:,1:4);

else
    error('mshReadMesh : unavailable case');
end

% Close file
fclose(fid);
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function BPI = myFilteredBackprojectionCentralSlice(sinogram,thetas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% filtered back projection in the frequency domain applying central slice
% theorem-> schlegel & bille 9.2.3
% modified by: Mark Bangert
% m.bangert@dkfz.de 2011
%
% note: a) matlab puts the 0 frequency component of a fourier spectrum _not_
%          in the middle. we need to fumble around with fftshift.

% figure out how big our picture is going to be.
numOfParallelProjections = size(sinogram,1);
numOfAngularProjections  = length(thetas);

% convert thetas to radians
thetas = (pi/180)*thetas;

% set up the backprojected image
BPI = zeros(numOfParallelProjections,numOfParallelProjections);

% find the middle index of the projections
midindex = floor(numOfParallelProjections/2) + 1;

% set up the coords of the image
[xCoords,yCoords] = meshgrid(ceil(-numOfParallelProjections/2):ceil(numOfParallelProjections/2-1));

% set up filter
rampFilter = [floor(numOfParallelProjections/2):-1:0 1:ceil(numOfParallelProjections/2-1)]';

% loop over each projection
for i = 1:numOfAngularProjections

    % figure out which projections to add to which spots
    rotCoords = round(midindex + xCoords*sin(thetas(i)) + yCoords*cos(thetas(i)));

    % check which coords are in bounds
    indices   = find((rotCoords > 0) & (rotCoords <= numOfParallelProjections));
    newCoords = rotCoords(indices);

    % filter
    filteredProfile = real( ifft( ifftshift( rampFilter.*fftshift(fft(sinogram(:,i)) ) ) ) );

    % summation
    BPI(indices) = BPI(indices) + filteredProfile(newCoords)./numOfAngularProjections;

    % visualization on the fly
    imagesc(BPI)
    drawnow

end


text:
%% Regression test data generator
% This program should be executed at the very beginning of refactoring a
% major revision. This is for making legacy STRAIGHT v40_007d to be
% compatible with MATLAB R2015b and Octave

% Copyright(c) 2016, Hideki Kawahara, (kawahara@sys.wakayama-u.ac.jp)

clear all
close all

original_speech_dir = '~/Music/VCTK_CORPUS/VCTK-Corpus/wav48/';
target_analysis_dir = '~/m-file/STRAIGHTV40_007e/analysisData/';
target_wave_dir = '~/m-file/STRAIGHTV40_007e/waveData/';

mkdir(target_analysis_dir);
mkdir(target_wave_dir);

dir_list = dir([original_speech_dir 'p*']);

%%

n_dirs = length(dir_list);
n_files = 0;
for ii = 1:n_dirs
  tmp_files = dir([original_speech_dir dir_list(ii).name '/*.wav']);
  n_files = n_files + length(tmp_files);
end;
%%
n_test = 2; % number of files tested for each speaker
l_segment = 0.1; % 100 ms segment
for ii = 1:n_dirs
  rng(12345); % initialize frozen random number
  seq_id = 0;
  basic_stat_table = zeros(n_files, 4);
  tmp_files = dir([original_speech_dir dir_list(ii).name '/*.wav']);
  for jj = 1:length(tmp_files)
    [x, fs] = audioread([original_speech_dir dir_list(ii).name '/' ...
      tmp_files(jj).name ]);
    seq_id = seq_id + 1;
    l_in_sample_segment = min(length(x), round(fs * l_segment));
    n_segment = floor(length(x) / l_in_sample_segment);
    rms_level = zeros(n_segment, 1);
    for kk = 1:n_segment
      rms_level(kk) = 20 * ...
        log10(std(x((kk - 1) * l_in_sample_segment + ...
        (1:l_in_sample_segment))));
    end;
    sorted_level = sort(rms_level);
    basic_stat_table(seq_id, 1) = length(x) / fs;
    basic_stat_table(seq_id, 2) = max(rms_level) - min(rms_level);
    basic_stat_table(seq_id, 3) = ...
      sorted_level(round(length(sorted_level) * 0.85));
    basic_stat_table(seq_id, 4) = max(abs(x));
  end;
  basic_stat_table = basic_stat_table(1:seq_id, :);
  % select safe region
  sorted_length = sort(basic_stat_table(:, 1));
  sorted_dynamic_range = sort(basic_stat_table(:, 2));
  sorted_85percent = sort(basic_stat_table(:, 3));
  l_10 = sorted_length(round(seq_id * 0.1));
  l_90 = sorted_length(round(seq_id * 0.9));
  d_10 = sorted_dynamic_range(round(seq_id * 0.1));
  d_90 = sorted_dynamic_range(round(seq_id * 0.9));
  v_10 = sorted_85percent(round(seq_id * 0.1));
  v_90 = sorted_85percent(round(seq_id * 0.9));
  index_list = 1:seq_id;
  safe_index = index_list( ...
    l_10 < basic_stat_table(:, 1) & l_90 > basic_stat_table(:, 1) & ...
    d_10 < basic_stat_table(:, 2) & d_90 > basic_stat_table(:, 2) & ...
    v_10 < basic_stat_table(:, 3) & v_90 > basic_stat_table(:, 3) & ...
    basic_stat_table(:, 4) < 0.95);
  selection_index = 1:length(safe_index);
  [~, tmp_index] = sort(rand(n_test, 1));
  selection_index = selection_index(tmp_index(1:n_test));
  for kk = 1:n_test
    id = safe_index(selection_index(kk));
    [x, fs] = audioread([original_speech_dir dir_list(ii).name '/' ...
      tmp_files(id).name ]);
    rng(12345); % initialize frozen random number
    f0raw = MulticueF0v14(x,fs);
    ap = exstraightAPind(x,fs,f0raw);
    n3sgram=exstraightspec(x,f0raw,fs);
    rng(12345); % initialize frozen random number
    y = exstraightsynth(f0raw,n3sgram,ap,fs);
    disp([num2str(kk) ': ' tmp_files(id).name  ' at:' datestr(now)]);
    audiowrite([target_wave_dir '/' tmp_files(id).name], ...
      y / max(abs(y)) * 0.9, fs);
    path_name_f0 = [target_analysis_dir '/' tmp_files(id).name(1:end-4) 'f0.bin'];
    path_name_ap = [target_analysis_dir '/' tmp_files(id).name(1:end-4) 'ap.bin'];
    path_name_sp = [target_analysis_dir '/' tmp_files(id).name(1:end-4) 'sp.bin'];
    WriteBinaryData(path_name_f0, f0raw)
    WriteBinaryData(path_name_ap, ap)
    WriteBinaryData(path_name_sp, n3sgram)
  end;
end;

author: HidekiKawahara


text:
function dat = read_spmeeg_data(filename, varargin)

% read_spmeeg_data() - import SPM5 and SPM8 meeg datasets
%
% Usage:
%   >> header = read_spmeeg_data(filename, varargin);
%
% Inputs:
%   filename - [string] file name
%
% Optional inputs:
%   'begsample'      first sample to read
%   'endsample'      last sample to read
%   'chanindx'  -    list with channel indices to read
%   'header'    - FILEIO structure header
%
% Outputs:
%   dat    - data over the specified range
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Vladimir Litvak

if nargin < 1
    help read_spmeeg_data;
    return;
end;

typenames = {'uint8','int16','int32','float32','float64','int8','uint16','uint32'};
typesizes   = [1  2  4  4 8 1 2 4];

header    = ft_getopt(varargin, 'header');
begsample = ft_getopt(varargin, 'begsample');
endsample = ft_getopt(varargin, 'endsample');
chanindx  = ft_getopt(varargin, 'chanindx');

if isempty(header)
    header = read_spmeeg_header([filename(1:(end-3)) 'mat']);
end

if isempty(begsample), begsample = 1; end;
if isempty(endsample), endsample = header.nSamples; end;

datatype = 'float32-le';
scale = [];
if isfield(header, 'orig')
    if isfield(header.orig, 'data') && isnumeric(header.orig.data) ...
            && ~isempty(header.orig.data)
        try
            dat = reshape(header.orig.data(chanindx, :, :), length(chanindx), []);
            dat = dat(:, begsample:endsample);
            return;
        end
    end

    if isfield(header.orig, 'datatype')
        datatype = header.orig.datatype;
    elseif isfield(header.orig.data, 'datatype')
        datatype = header.orig.data.datatype;
    end
    if isfield(header.orig, 'scale')
        if isnumeric(header.orig.scale)
            scale = header.orig.scale;
        else
            scale = header.orig.scale.values;
        end
    elseif isfield(header.orig.data, 'scale')
        scale = header.orig.data.scale;
    end
end

stepsize = typesizes(strmatch(strtok(datatype, '-'), typenames));

filename = [filename(1:(end-3)) 'dat'];

fid = fopen_or_error(filename, 'r');
fseek(fid, stepsize*header.nChans*(begsample-1), 'bof');
[dat, siz] = fread(fid, [header.nChans, (endsample-begsample+1)], strtok(datatype, '-'));
fclose(fid);

if ~isempty(chanindx)
    % select the desired channels
    dat = dat(chanindx,:);
end

if ~isempty(scale) && ~ismember(strtok(datatype, '-'), {'float32', 'float64'})

    % This is a somewhat complicated mechanism to figure out which scaling
    % coefficients go with which data points in a generic way

    trlind = floor(((begsample:endsample)-1)/header.nSamples)+1;

    utrlind = unique(trlind);

    for i = 1:length(utrlind)
        dat(:, trlind == utrlind(i)) = dat(:, trlind == utrlind(i)).* ...
            repmat(squeeze(scale(chanindx,utrlind(i))), 1, sum(trlind == utrlind(i)));
    end

end

author: spm


text:
function plotIntCon(prob,data)
%PLOTINTCON Plot Integer Constraints on the current figure
%   plotIntCon(prob)

%   Copyright (C) 2011 Jonathan Currie (I2C2)

%Check we have integer variables in this plot
if(all(prob.int.str(data.idx) == 'C'))
    return; %nothing to do!
end

xl = round(xlim); yl = round(ylim);
hold on;

%Get Problem Components
if(~isempty(prob.rl))
    [A,b,Aeq,beq] = row2gen(prob.A,prob.rl,prob.ru);
else
    A = prob.A; b = prob.b;
    Aeq = prob.Aeq; beq = prob.beq;
end
Q = prob.Q; l = prob.l; qrl = prob.qrl; qru = prob.qru;
if(isempty(prob.lb)), lb = -Inf; else lb = prob.lb; end
if(isempty(prob.ub)), ub = Inf; else ub = prob.ub; end

%Check for 1D problem
if(prob.sizes.ndec == 1)
    if(prob.int.ind(1) == -1)
        X = 0:1; Y = [];
    else
        X = (floor(xl(1))+1):(ceil(xl(2))-1); Y = [];
    end
else
    %Generate Grid and check for binary constraints
    if(prob.int.ind(1) == -1)
        x = 0:1;
    else
        try
            x = (floor(xl(1))+1):(ceil(xl(2))-1);
        catch ME
            optiwarn('opti:plotint',['Could not plot integer constraints due to plotting error: ' ME.message]);
            return;
        end
    end
    if(prob.int.ind(2) == -1)
        y = 0:1;
    else
        try
            y = (floor(yl(1))+1):(ceil(yl(2))-1);
        catch ME
            optiwarn('opti:plotint',['Could not plot integer constraints due to plotting error: ' ME.message]);
            return;
        end
    end
    [X,Y] = meshgrid(x,y);
end

%Check for row nonlinear constraints
if(~isempty(prob.cl))
    prob = nrow2mix(prob,0,false);
end

%Work out feasible integer points
idx = zeros(numel(X),1); Xv = X(:); Yv = Y(:);
for i=1:numel(X)
    if(isempty(Yv))
        xy = Xv(i);
    else
        xy = data.fixval;
        xy(data.idx) = [Xv(i), Yv(i)]';
    end
    if(~isempty(prob.A))
        idx(i) = all(b >= A*xy) && all(xy >= lb) && all(xy <= ub);
    else
        idx(i) = all(xy >= lb) && all(xy <= ub);
    end
    if(~isempty(Aeq))
        idx(i) = idx(i) && (abs(Aeq*xy-beq) < 1e-6);
    end
    if(~isempty(Q))
        for n = 1:prob.sizes.nqc
            if(iscell(Q))
                nQ = Q{n}; nl = l(:,n); nrl = qrl(n); nru = qru(n);
            else
                nQ = Q; nl = l; nrl = qrl; nru = qru;
            end
            if(~isinf(nru))
                qu = all(xy'*nQ*xy + nl'*xy <= nru);
            else
                qu = true;
            end
            if(~isinf(nrl))
                ql = all(xy'*-nQ*xy - nl'*xy <= -nrl);
            else
                ql = true;
            end
            idx(i) = idx(i) && ql && qu;
        end
    end
    if(~isempty(prob.nle))
        vals = prob.nlcon(xy);
        leI = find(prob.nle == -1);
        geI = find(prob.nle == 1);
        eqI = find(prob.nle == 0);
        nlrhs = prob.nlrhs;
        idx(i) = idx(i) && all(vals(leI) <= nlrhs(leI)) && all(vals(geI) >= nlrhs(geI)) && all(abs(vals(eqI)-nlrhs(eqI)) < 1e-6);
    end
end
%Plot points
idx = logical(idx);
if(isempty(Y))
    Y = zeros(size(X));
    for i = 1:length(Y)
        Y(i) = prob.objective(X(i));
    end
end
plot(X(idx),Y(idx),'bo','markerface','b','markersize',5);
plot(X(~idx),Y(~idx),'bo','markersize',5);

hold off;

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef Source < handle & matlab.mixin.Heterogeneous
	% Source is a superclass for all electric current sources J.

	properties (SetAccess = immutable)
		lgrid  % {x_array, y_array, z_array}: locations (intercepts with normal axis) of grid planes of grid type "this.gt"
		laltgrid  % {x_array, y_array, z_array}: locations (intercepts with normal axis) of grid planes of grid type "alter(this.gt)"
		shape  % shape of source; used to draw source
		forceprim  % true or false: flag to control the behavior of get.l()s
	end

	properties (SetAccess = private)
		gt  % GT.prim or GT.dual: type of grid line along which dipole sources of this source are defined
	end

	properties (Dependent, SetAccess = immutable)
		l
	end

	methods (Abstract = true)
		[index_cell, Jw_patch] = generate_kernel(this, w_axis, grid3d)
	end

	methods
		function this = Source(lgrid_cell, laltgrid_cell, shape, forceprim)
			chkarg(istypesizeof(lgrid_cell, 'realcell', [1 Axis.count], [1 0]), ...
				'"lgrid_cell" should be length-%d row cell array whose each element is row vector with real elements.', Axis.count);
			chkarg(istypesizeof(laltgrid_cell, 'realcell', [1 Axis.count], [1 0]), ...
				'"laltgrid_cell" should be length-%d row cell array whose each element is row vector with real elements.', Axis.count);
			chkarg(istypesizeof(shape, 'Shape'), '"shape" should be instance of Shape.');

			if nargin < 4  % no forceprim
				forceprim = false;
			end
			chkarg(istypesizeof(forceprim, 'logical'), '"forceprim" should be logical.');

			this.lgrid = lgrid_cell;
			this.laltgrid = laltgrid_cell;
			this.shape = shape;
			this.gt = GT.empty();
			this.forceprim = forceprim;
		end

		function set_gridtype(this, gt)
			% If this is SRCJ and E-field grid is primary, then gt = GT.prim.
			% If this is SRCJ and E-field grid is dual, then gt = GT.dual.
			% If this is SRCM and E-field grid is primary, then gt = GT.dual.
			% If this is SRCM and E-field grid is dual, then gt = GT.prim.
			chkarg(istypesizeof(gt, 'GT'), '"gt" should be instance of GT.');
			this.gt = gt;
		end

		function l = get.l(this)
			l = cell(Axis.count, GT.count);
			if this.forceprim
				l(:, GT.prim) = this.lgrid.';
				l(:, GT.dual) = this.laltgrid.';
			else
				l(:, this.gt) = this.lgrid.';
				l(:, alter(this.gt)) = this.laltgrid.';
			end
		end

		function [index_cell, JMw_patch] = generate(this, w_axis, grid3d)
			chkarg(istypesizeof(w_axis, 'Axis'), '"w_axis" should be instance of Axis.');
			chkarg(istypesizeof(grid3d, 'Grid3d'), '"grid3d" should be instance of Grid3d.');

			try
				[index_cell, JMw_patch] = this.generate_kernel(w_axis, grid3d);  % Cw_patch: current source
			catch err
				exception = MException('Maxwell:srcAssign', 'Source assignment failed.');
				throw(addCause(exception, err));
			end
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function y = randsample(n,k,j)
%function y = randsample(n,k)
%
%sample without replacement k integers from 1:n
if nargin<2, k = n; end
if nargin<3, j = 1; end
if numel(n) == 1
    n = 1:n;
end
[ig idx] = sort(rand(numel(n),j));
y = n(idx(1:k,:));
end

author: jacobeisenstein

repo: SAGE


text:
function P=rayleighpdf(M,S)

P= (((M./(S.^2)).*exp((-1.*(M.^2))./(2.*S^2)))).*(M>=0);

%%
% _*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 e = prtRvUtilDiscreteEntropy(q)
% DISCRETEENTROPY

e = q.*log(q);
e(q==0) = 0;
e = sum(e(:));

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function f = logical(f)
%LOGICAL   CLASSICFUN logical.
%   LOGICAL(F) returns a CLASSICFUN which evaluates to one at all points where F is
%   nonzero and zero otherwise.  F cannot have any roots in its domain.  If F
%   does have roots, then LOGICAL(F) will return garbage with no warning. F may
%   be complex.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

f.onefun = logical(f.onefun);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [rxnList,metList] = autoMatchReactions(scoreTotal, rxnList, metList, varargin)
% Finds confident matches from CMODEL to TMODEL.
% reactions which can not be confidently assigned are reviewed with
% `reactionCompare`. To make a confident assignment a matched reaction must
% have a score above a supplied cutoff and not have competing matches with
% similar scores. Called by `reactionCompareGUI`, calls `metCompare`.
%
% USAGE:
%
%    rxnList = autoMatchReactions(scoreTotal, rxnList, metList, [rxnHighCutoff, rxnMargin, rxnLowCutoff], [metHighCutoff, metMargin, metLowCutoff])
%
% INPUTS:
%    scoreTotal:       Optimized score matrix between `CMODEL` and `TMODEL` reactions.
%    rxnList:          Only undeclared reactions will be automatched.
%    metList:          Comparison array for metabolites
%
% OPTIONAL INPUTS:
%    rxnHighCutoff:    Score value above which confident assignments can be made.
%    rxnMargin:        Distance a best match has to be above a competitor match.
%    rxnLowCutoff:     Score value below which reactions will be declared new.
%    metHighCutoff:    The following three are the same but for mets.
%    metMargin:
%    metLowCutoff:
%    CMODEL:           global input
%    TMODEL:           global input
%
% OUTPUTS:
%    rxnList:          Array which correlates a reaction from CMODEL to the index
%                      of the best match in `TMODEL`. Rxns with 0 designation are
%                      new, ones with a -1 need manual review.
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

h = waitbar(0, 'Automatching') ; % Declare variables.

global CMODEL TMODEL

% Number of reactions in CMODEL.
nCrxn = size(scoreTotal, 1) ;

% Assign cutoff params if possible.
if nargin >= 6 %8
    rxnHighCutoff = varargin{1} ;
    rxnMargin = varargin{2} ;
    rxnLowCutoff = varargin{3} ;
else
    rxnHighCutoff = 0.99 ;
    rxnMargin = 0.1 ;
    rxnLowCutoff = 0.01 ;
end
% get metabolite automatch cutoffs if possible
if nargin >= 8 %10
    metHighCutoff = varargin{4} ;
    metMargin = varargin{5} ;
    metLowCutoff = varargin{6} ;
else
    metHighCutoff = 0.99 ;
    metMargin = 0.1 ;
    metLowCutoff = 0.01 ;
end
% Save a form of the rxnList before auto matching.
rxnListBefore = rxnList ;

%% Find probable matches based on scoring.
% Find best scoring matches.
[sortMatch, sortMatchIndex] = sort(scoreTotal, 2, 'descend') ;

% Only consider reactions which have not been previously declared.
probMatch = [] ;
pos = 1 ;
for iRxn = 1:nCrxn
    if sortMatch(iRxn, 1) > rxnHighCutoff && ...
            sortMatch(iRxn, 1) - sortMatch(iRxn, 2) > rxnMargin && ...
            rxnList(iRxn) == -1
        rxnList(iRxn) = sortMatchIndex(iRxn, 1) ;
        probMatch(pos, 1) = iRxn ;
        pos = pos + 1 ;
    end
end

%% Double check reactions for hints that they are misses.
% Set the test columns in probMatch to 1. Innocent until proven guilty.
if ~isempty(probMatch)
    probMatch(:, 2:5) = 1 ;
end

for iRxn = 1:size(probMatch, 1)
    % Check if matched reaction has been associated with another reaction.
    nowMatch = rxnList(probMatch(iRxn, 1)) ;
    % The match should only be in rxnList once.
    if length(find(rxnList == nowMatch)) ~= 1
        probMatch(iRxn, 2) = 0 ;
    end

    % Check if paired reactions have KEGG IDs and if they match.
    cKegg = CMODEL.rxnKEGGID{probMatch(iRxn, 1)} ;
    tKegg = TMODEL.rxnKEGGID{rxnList(probMatch(iRxn, 1))} ;
    if ~isempty(cKegg) && ~isempty(tKegg)
        pipePos = [0 strfind(cKegg, '|') length(cKegg) + 1] ;
        for jID = 1:length(pipePos) - 1
            nowID = cKegg(pipePos(jID)+1:pipePos(jID + 1) - 1) ;
            if isempty(strfind(tKegg, nowID))
                probMatch(iRxn, 3) = 0 ;
            end
        end
    end

    % Check if paired reactions have the same stoichiometry.
    if CMODEL.metNums(probMatch(iRxn, 1), 2) ~= ...
            TMODEL.metNums(rxnList(probMatch(iRxn, 1)), 2)
        probMatch(iRxn, 4) = 0 ;
    end

    % Check if paired reactions have the same compartment.
    if ~isempty(CMODEL.rxnComp{probMatch(iRxn, 1)} ) && ~isempty(TMODEL.rxnComp{rxnList(probMatch(iRxn, 1))})
        if ~strcmp(CMODEL.rxnComp{probMatch(iRxn, 1)}, ...
                TMODEL.rxnComp{rxnList(probMatch(iRxn, 1))})
            probMatch(iRxn, 5) = 0 ;
        end
    end

    % If any of the above tests failed, set the match in rxnList backto -1.
    if probMatch(iRxn, 2:end)
    else
        rxnList(probMatch(iRxn, 1)) = -1 ;
    end
end

clear probMatch nowMatch cKegg tKegg

%% Assign new reactions based on cutoff.
for iRxn = 1:nCrxn
    if rxnList(iRxn) == -1 && ...
            sortMatch(iRxn,1) < rxnLowCutoff
        rxnList(iRxn) = 0 ;
    end
end

%% Assign new reactions based on metabolites that have been declared.
% For undeclared reactions in CMODEL, find the indexes of their metabolite
% matches.
for iRxn = 1:length(rxnList)
    if rxnList(iRxn) == -1
        % For the reaction in question from CMODEL.
        metPos = find(CMODEL.S(:, iRxn)) ;
        cMetMatchIndexes = zeros(length(metPos), 1) ;
        for iMet = 1:length(metPos)
            cMetMatchIndexes(iMet) = metList(metPos(iMet)) ;
        end
        % For the best match from TMODEL.
        tMetMatchIndexes = find(TMODEL.S(:, sortMatchIndex(iRxn, 1))) ;
        % If all the mets have been declared in cRxn, and the best match
        % from T does not have exactly those mets, declare a new reaction.
        if isempty(find(cMetMatchIndexes == 0, 1))
            if logical(ismember(cMetMatchIndexes, tMetMatchIndexes)) * ...
                    length(cMetMatchIndexes) == length(tMetMatchIndexes)
            else
                rxnList(iRxn) = 0 ;
            end
        end
        % Also declare reactions as new it they have a new metabolite(s).
        if ~isempty(find(cMetMatchIndexes > length(TMODEL.mets), 1))
            rxnList(iRxn) = 0 ;
        end
    end
end

%% For reactions that have just been matched or declared new, find mets.
% Find newly auto-matched reactions.
autoMatchedRxns = find(rxnList(:) ~= rxnListBefore(:)) ;

% Prepare data structure to be passed.
RxnInfo.rxnList = rxnList ;

for iRxn = 1:length(autoMatchedRxns)
    RxnInfo.rxnIndex = autoMatchedRxns(iRxn) ;
    RxnInfo.rxnMatch = rxnList(autoMatchedRxns(iRxn)) ;
    RxnInfo.metList = metList ;
    RxnInfo.metAutoMatchLimits = [metHighCutoff, metMargin, metLowCutoff] ;

    % Launch comparison script.
    [metList, stopFlag] = metCompare(RxnInfo) ;

    % If metCompare was suspended, don't attempt to find matches for the
    % reamining reactions.
    if stopFlag
        break
    end
    waitbar(iRxn/length(autoMatchedRxns), h)
end

close(h)

author: opencobra

repo: cobratoolbox


text:
function [Y,W,SetupStruc] = Process_MVDR_Search(s,Transfer,SetupStruc)
K = SetupStruc.MVDR_Search.K;
hop = SetupStruc.MVDR_Search.hop;
win = hanning(K,'periodic');
win = win/sqrt(sum(win(1:hop:K).^2));
SetupStruc.MVDR_Search.win = win;  % Preserve 'win' in 'SetupStruc'
% K = SetupStruc.MVDR.K;
% hop = SetupStruc.MVDR.hop;
% win = hanning(K,'periodic');
% win = win/sqrt(sum(win(1:hop:K).^2));
% SetupStruc.MVDR.win = win;  % Preserve 'win' in 'SetupStruc'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(s,2);
for i = 1:N
    X(:,:,i) = fft(enframe(s(:,i),win,hop)');
end
frame_N = size(X,2);
K_m = K/2+1;
Num = size(Transfer,3);
Y = zeros((frame_N-1)*hop+K,Num);
Y_f = zeros(size(X,1),size(X,2),Num);
WNG = zeros(K_m,Num);
%%%%%%%%%%%%%%%%%%%%%%%%%% Obtain processing matrix 'W'
theta = 10^-10;
WNG_set = 0;
W = zeros(Num,N,K_m);
for i = 2:K_m
    X_f = permute(X(i,:,:),[3 2 1]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% MVDR processing
    W_f = zeros(N,Num);
    Steer = permute(Transfer(i,:,:),[2 3 1]);
    for j = 1:Num
        R = X_f*X_f'/frame_N;
        if rcond(R)<theta
            R = R+eye(N)*min([min(diag(R)) theta]);
        end
        h = Steer(:,j);
        w = R\h/(h'/R*h);
        wng = -10*log10(diag(w'*w));
        e = 0.05;
        if(WNG_set-wng>0)
            sign = 1;
        else
            sign = 0;
        end
        signM = 1;
        while(sign==1 && abs(wng-WNG_set)>0.5)
            if(wng>WNG_set)
                R = R-e*eye(N);
                e = e/10;
                signM = 0;
            else
                if(signM==1)
                    e = e*2;
                end
                R = R+e*eye(N);
            end
            w = R\h/(h'/R*h);
            wng = -10*log10(diag(w'*w));
        end
%         while(sign==0 && abs(wng-WNG_set)>0.5)
%             if(wng<WNG_set)
%                 R = R+e*eye(N);
%                 e = e/10;
%                 signM = 0;
%             else
%                 if(signM==1)
%                     e = e*2;
%                 end
%                 R = R-e*eye(N);
%             end
%             w = R\h/(h'/R*h);
%             wng = -10*log10(diag(w'*w));
%         end
        W_f(:,j) = w;
    end
    W_f = W_f';
    W(:,:,i) = W_f;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    WNG(i,:) = -10*log10(diag(W_f*W_f')');
    Y_ = W_f*X_f;
    Y_f(i,:,:) = Y_.';
    if(i~=K_m)
        Y_f(K+2-i,:,:) = Y_';
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recover signals
if(K/hop==2)
    win = ones(K,1);
end
for i = 1:Num
    Y(:,i) = overlapadd(real(ifft(Y_f(:,:,i)))',win,hop);
end
% autoPlot(WNG,'MVDR_Search');
return;

repo: Voice-Separation-and-Enhancement


text:
function [y,deriv] = scaleColumns(w,scales)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
%   w --> bsxfun(@times,reshape(w,[],n),scales(:)')
%
%   where n = length(scales);
%
% Note: this is a symmetric linear transform.

if nargin==0
    test_this();
    return;
end

if isempty(w)
    map = @(w)map_this(w,scales);
    y = linTrans(w,map,map);
    return;
end

if isa(w,'function_handle')
    f = scaleColumns([],scales);
    y = compose_mv(f,w,[]);
    return;
end

f = scaleColumns([],scales);
if nargout==1
    y = f(w);
else
    [y,deriv] = f(w);
end

function w = map_this(w,scales)
n = length(scales);
w = reshape(w,[],n);
w = bsxfun(@times,w,scales(:)');

function test_this()
K = 5;
N = 10;
M = randn(K,N);
scales = randn(1,N);

f = scaleColumns([],scales);
test_MV2DF(f,M(:));

author: nesl


text:
function Res = AccSumK(p,K)
%ACCSUMK      Res represents K-fold faithful rounding of sum(p)
%
%   Res = AccSumK(p,K)
%
%On return, Res represents a K-fold faithful rounding of sum(p), also in the
%  presence of underflow. Input vector p may be single or double precision.
%
%Implements Algorithm 6.4 from
%  S.M. Rump, T. Ogita, S. Oishi: Accurate Floating-point Summation II:
%    Sign, K-fold Faithful and Rounding to Nearest, Siam J. Sci. Comput.,
%    31(2):1269-1302, 2008.
%Requires (4m+5K+3)n flops for m executions of repeat-until loop in the
%  first and one execution of the repeat-until loops in subsequent calls
%  of TransformK.
%
%Reference implementation! Slow due to interpretation!
%

% written  03/03/07     S.M. Rump
% modified 05/09/09     S.M. Rump  rounding to nearest, complex input
%

  if ~isreal(p)
    Res = complex(AccSumK(real(p),K),AccSumK(imag(p),K));
    return
  end

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

  if isa(p,'double')
    nmax = 2^26-2;          % nmax = 67,108,864
  else
    nmax = 2^12-2;          % nmax = 4,094
  end
  if length(p)>nmax
    error(['maximum length of input vector for AccSumK ' int2str(nmax) '.'])
  end

  R = 0;
  if isa(p,'double'), prec='double'; else prec='single'; end
  Res = zeros(1,K,prec);
  [Res(1),R,p,sigma,Ms] = TransformK(p,R);
  if abs(Res(1))<=realmin(prec)
    if rndold, setround(rndold); end
    return
  end
  for k=2:K
    [Res(k),R,p,sigma,Ms] = TransformK(p,R,sigma,Ms);
    if abs(Res(k))<=realmin(prec)
      if rndold, setround(rndold); end
      return
    end
  end

  if rndold
    setround(rndold)
  end


text:
close all;
num = size(AllFeature1,2);
F1 = AllFeature1' / 4000;
% F1 = max(AllFeature1(1:512,:)', AllFeature1(513:1024,:)');%(1:512,:)
% F1 = bsxfun(@minus, F1, mean(F1,1));
% F1 = bsxfun(@rdivide, F1, sqrt(sum(F1.^2,2)));

F2 = AllFeature2' / 4000;
% F2 = max(AllFeature2(1:512,:)', AllFeature2(513:1024,:)');%(1:512,:)
% F2 = bsxfun(@minus, F2, mean(F2,1));
% F2 = bsxfun(@rdivide, F2, sqrt(sum(F2.^2,2)));

thresh2 = zeros(num,1);
for i = 1:num
%     thresh2(i) = F1(i,:) * mapping.A * F1(i,:)' + F2(i,:) * mapping.A * F2(i,:)' - 2 * F1(i,:) * mapping.G * F2(i,:)';
    thresh2(i) = pdist2(F1(i,:) ./ norm(F1(i,:)),F2(i,:) ./ norm(F2(i,:)));
%     thresh2(i) = F1(i,:) * F2(i,:)';
end;
figure;
hist(thresh2(1:3000),500);
figure;
hist(thresh2(3001:end),500);

accuracies = zeros(10,1);
for i=1:10
    test_idx = [(i-1) * 300 + 1 : i*300, (i-1) * 300 + 3001 : i*300 + 3000];
    train_idx = 1:6000;
    train_idx(test_idx) = [];
    bestc=256;
    same_label = ones(6000,1);
    same_label(3001:6000) = 0;

    mean_feature = mean([F1(train_idx,:); F2(train_idx,:)]);
    std_feature = std([F1(train_idx,:); F2(train_idx,:)]);
    F1_mu = bsxfun(@minus, F1, mean_feature);
    F2_mu = bsxfun(@minus, F2, mean_feature);
%     F1_mu = bsxfun(@rdivide, F1_mu, std_feature);
%     F2_mu = bsxfun(@rdivide, F2_mu, std_feature);
    F1_mu = bsxfun(@rdivide, F1_mu, sqrt(sum(F1_mu.^2,2)));
    F2_mu = bsxfun(@rdivide, F2_mu, sqrt(sum(F2_mu.^2,2)));
    F1_mu = F1;
    F2_mu = F2;

    [U,mu,vars] = pca( [F1_mu(train_idx,:); F2_mu(train_idx,:)]' );
    sum_var = cumsum(vars);
    sum_var = sum_var / sum_var(end);
    dims = find(sum_var > 0.995, 1, 'first');
    [F1PCA,F1PCAHat,~] = pcaApply( F1_mu', U, mu, dims );
    [F2PCA,F2PCAHat,~] = pcaApply( F2_mu', U, mu, dims );
    F1PCA = bsxfun(@rdivide, F1PCA, sqrt(sum(F1PCA.^2)));
    F2PCA = bsxfun(@rdivide, F2PCA, sqrt(sum(F2PCA.^2)));

    thresh1 = zeros(num,1);
    for n = 1:num
%         thresh1(n) = pdist2(F1_mu(n,:),F2_mu(n,:));
%         thresh1(n) = pdist2(F1(n,:),F2(n,:));
        thresh1(n) = pdist2(F1PCA(:,n)',F2PCA(:,n)');
%         thresh1(n) = F1PCA(:,n)' * F2PCA(:,n);
%         thresh1(n) = F1(n,:) * F2(n,:)';
    end;
    cmd = [' -t 0 -h 0'];
    model = svmtrain(same_label(train_idx),thresh1(train_idx),cmd);
    [class, accuracy, deci] = svmpredict(same_label(test_idx),thresh1(test_idx),model);
%     [~, threshold_fold] = Sys_accuracy(thresh1(train_idx(1:2700)), thresh1(train_idx(2701:5400)));
%     accuracy = (sum(thresh1(test_idx(1:300)) < threshold_fold) + sum(thresh1(test_idx(301:600)) >= threshold_fold)) / 600;
    fprintf('%d th-fold accuracy:%.4f\n', i, accuracy(1));
    accuracies(i) = accuracy(1);
end;
fprintf('(PCA)accuracy by 10-fold evalution:%.4f\r\n', mean(accuracies));
cmd = [' -t 0 -h 0'];
model = svmtrain(same_label,thresh2,cmd);
[class, accuracy, deci] = svmpredict(same_label,thresh2,model);
% mean(thresh2(same_label==1)) / 4 + mean(max(0,1 - thresh2(same_label==0))) / 4
% sum((thresh2<0.22) == same_label) / 6000
% [accuracy, threshold] = Sys_accuracy(thresh2(same_label==1), thresh2(same_label==0));
% fprintf('(Without PCA)accuracy by direct evalution:%.4f, threshold:%.4f\r\n', accuracy, threshold);

author: happynear

repo: FaceVerification


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
    = cornk(fid, data, varargins, opts)
% This program starts the CORN-K algorithm
% CORN-K is a top-K combination of CORN experts
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ....
%           = cornk_start(fid, data, varargins, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%          = cornk_start(fid, data, {5, 10, 0.1, 0}, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the parameters
K  = varargins{1};
L  = varargins{2};
pc  = varargins{3};      % percentage of top-k
tc = varargins{4};      % transaction cost fee rate

% Run the CORN-K algorithm
[cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
    = cornk_run(fid, data, K, L, pc, tc, opts);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function gIX = HierClus_Direct(C,gIX)
[gIX, numU] = SqueezeGroupIX(gIX);
% [C,~] = FindCentroid_Direct(gIX,M);
D = pdist(C,'correlation');
if size(C,1)>1,
    tree = linkage(C,'average','correlation');
    leafOrder = optimalleaforder(tree,D);

    % sort for uniform colorscale
    temp = zeros(size(gIX));
    for i = 1:numU,
        temp(gIX==leafOrder(i)) = i; % = T(i) for clusters segmented from tree
    end
    gIX = temp;
end
end

author: xiuyechen

repo: FishExplorer


text:
function [llh] = tapas_ti_llh(y, x, u, theta, ptheta)
%% General method to iterate across chains and subjects
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

fllh = ptheta.method_llh;

ns = size(theta, 1);
nc = size(theta, 2);

llh = zeros(1, nc);

for i = 1:nc
    tllh = zeros(ns, 1);
    for j = 1:ns
        tllh(j) = fllh(y{j, i}, x{j, i}, u{j, i}, theta{j, i}, ptheta);
    end
    llh(i) = sum(tllh);
end

%if any(llh < 1000)
%    keyboard
%end

end

author: translationalneuromodeling

repo: tapas


text:
function varargout = barvalues(h,precision)
% BARVALUES  Display bar values ontop of bars in bar or histogram plot.
% SYNTAX:
% barvalues;                  - operates on currnet axes.
% barvalues(h);               - operates on h.
% barvalues(_,precision);     - additionaly, specifies the precision of
%                               displayed values. or format.
% t = barvalues(_,_);         - returns the handles to the value text objects.
%
% h - handle to axes or bar (operates on specified object only)
%     or figure (operates on all child axes).
%
% precision - Decimal precision to display (0-10),
%             or 'formatSpec' as in num2str. (default:'% .0f')
%
% t - handles to the text objects.
%
%   For more information about 'formatSpec':
%   See also NUM2STR.

%Author: Elimelech Schreiber, 11/2017
% ver 2.1    - updated 04/2018

t=[];

if nargin>1 && ~isempty(precision) % Parse precision
    if isnumeric(precision) && precision >=0 && precision <=10
        precision =['% .',int2str(precision),'f'];
    elseif ~ischar(precision) && ~isstring(precision)
        error('Precision format unsupported.');
    end
else
    precision ='% .0f';
end

if nargin<1 || isempty(h)   % parse h (handle)
    h =gca;
elseif isaType(h,'figure')
   B = findobj(h,'type','bar','-or','type','Histogram'); % apply to multiple axes in figure.
   for b =B'
           t = [t; {barvalues(b,precision)}]; % Return array of text objects
                                              % for each bar plot.
   end
    if nargout>0
        varargout{1}=t;
    end
    return;
end
if isaType(h,'axes')
    h  = findobj(h,'type','bar','-or','type','Histogram','-or','type','patch');
    if isempty(h)
        return; % silently. to support multiple axes in figure.
    end
end
h = h(isaType(h,'bar') | isaType(h,'patch') | isaType(h,'histogram'));
if isempty(h)
    error('Cannot find bar plot.');
end
if size(h,1)>size(h,2)
    h=h';
end
for hn = h
    axes(ancestor(hn,'axes')); % make intended axes curent.
    if isaType(hn,'histogram')
        t =[t;  histvalues(hn,precision)];
        continue;
    end
    if isfield(hn,'XOffset')&&~isempty(hn.XOffset)
        XOffset = hn.XOffset;
    else
        XOffset = 0;
    end
    if isfield(hn,'YOffset')&&~isempty(hn.YOffset)
        YOffset = hn.YOffset;
    else
        YOffset = 0;
    end
    xData = hn.XData +XOffset;
    yData = hn.YData +YOffset;
    if size(xData,1)==4
      xData=mean(xData);
      yData=yData(2,:);
    end
    t = [t;  text(xData,yData,...    %position
       arrayfun(@(x)num2str(x,precision),yData,'UniformOutput' ,false),...    %text to display
        'HorizontalAlignment','center','VerticalAlignment','bottom')];
end
if nargout>0
    varargout{1}=t;
end

function flag =isaType(h,type)
try
    flag =strcmpi(get(h, 'type'), type);
catch
    flag =false;
end

function flag = isfield(h,fld)
flag =true;
try
    get(h,fld);
catch
    flag =false;
end

function t =histvalues(h,precision)
   hn=h;
    axes(ancestor(hn,'axes')); % make intended axes curent.
%     if isfield(hn,'XOffset')&&~isempty(hn.XOffset), XOffset = hn.XOffset; else XOffset = 0; end
%     if isfield(hn,'YOffset')&&~isempty(hn.YOffset), YOffset = hn.YOffset; else YOffset = 0; end
    xData = (hn.BinEdges(1:end-1) + hn.BinEdges(2:end))/2; yData = hn.Values;

    t = text(xData,yData,...    %position
       arrayfun(@(x)num2str(x,precision),yData,'UniformOutput' ,false),...    %text to display
        'HorizontalAlignment','center','VerticalAlignment','bottom');


text:
clear all;close all;clc

M_dir = GetFishDirectories();

varList = {'CR_dtr','nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn'};

i_fish = 10;

%% load data
disp(['load fish ' num2str(i_fish)]);
tic
datadir = M_dir{i_fish};
load(fullfile(datadir,['Fish' num2str(i_fish) '_direct_load_nodiscard.mat']),varList{:});
toc

CR_dtr_full = CR_dtr;
CInfo_full = CInfo;
%% Test: discard 50% noisy cells based on std of zscore of baseline

% Compute std of dimest 10% of frames for each cell.
tic
prc = prctile(CR_dtr_full,10,2);
nCells = size(CR_dtr_full,1);
STD_full = zeros(nCells,1);
for i = 1:nCells,
    ix = find(CR_dtr_full(i,:)<prc(i));
    STD_full(i) = std(CR_dtr_full(i,ix));
end
toc
beep

%% Set threshold at 50%, i.e. discard 50% of all cells
perc_keep = 50;
thr = prctile(STD_full,perc_keep);

temp = find(STD_full>thr);
[~,I] = sort(STD_full(temp));
cIX = temp(I);

%% visualize cells to discard
% gIX = (round((1:length(cIX))/1000)+1)'; % option: view sorted index
gIX = round(cIX/1000)+1; % option: view anatomical z index
M = CR_dtr_full(cIX,1:1000);
[~,stim] = StimulusKeyPreprocessing(frame_turn,i_fish);
BasicPlotMaps(cIX,gIX,M,CInfo,stim,anat_yx,anat_yz,anat_zx);

% I_v: index of valid cells
I_v_Holder = ones(1,nCells);
I_v_Holder(cIX) = 0;
I_v = find(I_v_Holder);

CR_dtr = CR_dtr_full(I_v,:);
STD = STD_full(I_v);
nCells = length(I_v);
CInfo = CInfo_full(I_v);

%% visualize valid cells, sorted by noise
[~,I] = sort(STD);
cIX = I;
gIX = (round((1:length(cIX))/1000)+1)'; % sorted index
M = CR_dtr(I,1:1000);
BasicPlotMaps(cIX,gIX,M,CInfo_thr,stim,anat_yx,anat_yz,anat_zx);

%% Save new 'directload.mat'
temp = fullfile(datadir,['Fish' num2str(i_fish) '_direct_load_thr50.mat']);
varList = {'CR_dtr','nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn','perc_keep'};
% NOTE! not saving CR_raw in non-'_nodiscard.mat' again!
% should keep '_nodiscard.mat' so that this thresholding can be re-done
save(temp,varList{:},'-v7.3');


text:
clear all; close all; clear classes; clc;

%% Set flags.
inspect_only = false;

%% Create shapes.
a = 420;  % lattice constant
t = 1;  % slab thickness
r = 0.29*a;  % hole radius

ad = 25;  % divider for a
td = 10;  % divider for t
dd = 10;  % divider for d = 2*r

mx = 20.5;
my = 7.5;
slab_yn = Box([-mx*a mx*a; -my*a -0.5*a; 0 t], [a/ad, a/ad, t]);
slab_yp = Box([-mx*a mx*a; 0.5*a my*a; 0 t], [a/ad, a/ad, t]);

rod = CircularCylinder(Axis.z, t, [0 0 t/2], r, [2*r/dd, 2*r/dd, t]);

%% Solve the system.
gray = [0.5 0.5 0.5];  % [r g b]
src_loc = 2*a;
[E, H, obj_array, src_array, J] = maxwell_run(...
	'OSC', 1e-9, 1550, ...
	'DOM', {'Johnson/Ag', gray}, [-mx*a mx*a; -my*a my*a; 0 t], [a/ad a/ad t], BC.p, [10*a 2*a 0], ...
	'OBJ', ...
		{'vacuum', 'w', 1}, Box([-mx*a mx*a; -a/2 a/2; 0 t]), ...
	'SRCJ', PointSrc(Axis.z, [src_loc, 0, 0.5]), ...
	inspect_only);

%% Visualize the solution.
if ~inspect_only
	figure;
	clear opts
% 	opts.withgrid = true;
	opts.withobjsrc = false;
	opts.withabs = true;
	opts.withpml = false;
	opts.phase = pi/2;
	vis2d(E{Axis.z}, Axis.z, 0.5, obj_array, src_array, opts)
	%%
	vis2d(H{Axis.y}, Axis.z, 0.5, obj_array, src_array, opts)

	%%
	flux_loc = 3*a;
	power_right = powerflux_patch(E, H, Axis.x, src_loc + flux_loc);
	power_left = -powerflux_patch(E, H, Axis.x, src_loc - flux_loc);
	fprintf('power:\n');
	fprintf('right = %s\n', num2str(power_right));
	fprintf('left = %s\n', num2str(power_left));
	fprintf('error = %s%%\n',num2str((power_left-power_right)/power_right*100));

	%%
	Sx = poynting(Axis.x, E{Axis.y}, E{Axis.z}, H{Axis.y}, H{Axis.z}, Axis.y, 0);
	[array, l] = Sx.data_original;
	plot(l{2}, abs(array))
	mx*a - 10*a
end


text:
function p = TimeExcludeRanges(tInc,t)

p=[];
mdiff = diff(tInc);
tidx1 = find(mdiff==-1);
tidx2 = find(mdiff==1);

% If excluded time points are at the extremes of the time range,
% then add inflections at the beginning or end.
if tInc(1)==0
    tidx1 = [1;tidx1];
end
if tInc(end)==0
    tidx2 = [tidx2;length(tInc)];
end
% if ~isempty(tidx2) && isempty(find(tidx1<tidx2(1)))
%     tidx1 = [1; tidx1];
% end
% if ~isempty(tidx1) && isempty(find(tidx2>tidx1(end)))
%     tidx2 = [tidx2; length(tInc)];
% end

for ii=1:length(tidx1)
    if ii<=length(tidx2)
        p(ii,:) = [t(tidx1(ii)) t(tidx2(ii))];
    else
        p(ii,:) = [t(tidx1(ii)) t(end)];
    end
end

author: BUNPC


text:
% This file create a simple regression based network for speech
% dereverberation or enhancement
%
% Created by Xiong Xiao, Temasek Laboratories, NTU, Singapore.
% Last Modified: 08 Feb 2017
%
function layer = genNetworkDereverb_Gaussian(para, stage)
para.freqBin = (0:1/para.fft_len:0.5)*2*pi; % w = 2*pi*f, where f is the normalized frequency k/N is from 0 to 0.5.
nFreqBin = length(para.freqBin);

% Part 1: generate the BF weight predicting subnet

if para.useWav  % input is row waveform
    layer{1}.name = 'Input';
    layer{end}.inputIdx = 1;
    layer{end}.dim = [1 1]*para.nCh;

    layer{end+1}.name = 'stft';
    layer{end}.prev = -length(layer)+1;
    layer{end}.fft_len = para.fft_len;
    layer{end}.frame_len = para.frame_len;
    layer{end}.frame_shift = para.frame_shift;
    layer{end}.removeDC = para.removeDC;
    layer{end}.win_type = para.win_type;
    layer{end}.dim = [(para.fft_len/2+1)*para.nCh layer{length(layer)+layer{end}.prev}.dim(1)];
    layer{end}.skipBP = 1;  % skip backpropagation

    layer{end+1}.name = 'Affine';       % scaling the Fourier transform
    layer{end}.prev = -1;
    layer{end}.W = [];
    layer{end}.b = [];
    layer{end}.dim = [1 1] * layer{length(layer)+layer{end}.prev}.dim(1);
    layer{end}.update = 0;
    layer{end}.skipBP = 1;  % skip backpropagation

    % extract a subset of dimensions for prediction. Not used currently, so we
    % will use all the dimensions.
    layer{end+1}.name = 'ExtractDims';
    layer{end}.prev = -1;
    layer{end}.dimIndex = 1:nFreqBin*para.nCh;
    layer{end}.dim = [length(layer{end}.dimIndex) layer{length(layer)+layer{end}.prev}.dim(1)];
    layer{end}.skipBP = 1;

    % get the log power spectrum and perform CMN
    layer{end+1}.name = 'Power';
    layer{end}.prev = -1;
    layer{end}.dim = [1 1] * layer{end-1}.dim(1);
    layer{end}.skipBP = 1;

    layer{end+1}.name = 'Log';
    layer{end}.const = 1e-2;
    layer{end}.prev = -1;
    layer{end}.dim = [1 1]*layer{end-1}.dim(1);
    layer{end}.skipBP = 1;
else    % input is log spectrogram
    layer{1}.name = 'Input';
    layer{end}.inputIdx = 1;
    layer{end}.dim = [1 1]*(para.fft_len/2+1)*para.nCh;
end

if para.useCMN
    layer{end+1}.name = 'CMN';
    layer{end}.prev = -1;
    layer{end}.dim = [1 1]*layer{end-1}.dim(1);
    layer{end}.skipBP = 1;
end

layer{end+1}.name = 'Delta';
layer{end}.prev = -1;
layer{end}.delta_order = 2;
layer{end}.dim = [layer{end}.delta_order+1 1]*layer{length(layer)+layer{end}.prev}.dim(1);
layer{end}.skipBP = 1;

layer{end+1}.name = 'Affine';       % scaling the Fourier transform
layer{end}.prev = -1;
layer{end}.W = [];
layer{end}.b = [];
layer{end}.dim = [1 1] * layer{length(layer)+layer{end}.prev}.dim(1);
layer{end}.update = 0;
layer{end}.skipBP= 1;

affineLayer = length(layer);

if ~isempty(para.hiddenLayerSizeShared)
    tmpTopology.inputDim = layer{end}.dim(1);
    tmpTopology.hiddenLayerSizeLSTM = para.hiddenLayerSizeShared;
    tmpTopology.usePastState = zeros(1,length(para.hiddenLayerSizeShared)); % do not use peeping hole
    tmpTopology.hiddenLayerSizeFF = [];
    tmpTopology.outputDim = 1024;
    tmpTopology.useAffineBtwLSTM = 1;
    tmpTopology.costFn = 'mse';
    tmpTopology.LastActivation4MSE = 'linear';
    layerShared = genNetworkLSTM(tmpTopology);
    if stage==2.2   % for stage 2.2, we don't even update the shared layer, only update the var subnet
        for i=1:length(layerShared)
            layerShared{i}.update = 0;
            layerShared{i}.skipBP = 1;
        end
    end
    layer = [layer layerShared(2:end-3)];
end

shared_lstm_idx = length(layer);

% get the mean
tmpTopology.inputDim = layer{shared_lstm_idx}.dim(1);
tmpTopology.hiddenLayerSizeLSTM = para.hiddenLayerSizeMu;
tmpTopology.usePastState = zeros(1,length(para.hiddenLayerSizeMu)); % do not use peeping hole
tmpTopology.hiddenLayerSizeFF = [];
if strcmpi(para.DeltaGenerationType, 'DeltaByEqn')
    tmpTopology.outputDim = nFreqBin;
else
    tmpTopology.outputDim = nFreqBin*3;
end
tmpTopology.costFn = 'mse';
tmpTopology.LastActivation4MSE = 'linear';
layerMu = genNetworkLSTM(tmpTopology);
layerMu = layerMu(2:end-2);
if ~isempty(para.hiddenLayerSizeMu) && ~isempty(para.hiddenLayerSizeShared)  % add a projection layer between two LSTM layers
    layerMu = [layer(affineLayer) layerMu];
    layerMu{1}.update = 1;
    layerMu{1}.skipBP = 0;
    layerMu{1}.dim = [1 1]*layer{shared_lstm_idx}.dim(1);
end
layerMu{1}.prev = shared_lstm_idx - length(layer)-1;
if floor(stage)==2
    for i=1:length(layerMu)
        layerMu{i}.update = 0;
        if stage==2.2  % if there is no shared layers, we don't let the gradient from mean layers
            layerMu{i}.skipBP = 1;
        end
    end
end
layer = [layer layerMu];

if strcmpi(para.DeltaGenerationType, 'DeltaByEqn')
    layer{end+1}.name = 'Delta';
    layer{end}.prev = -1;
    layer{end}.delta_order = 2;
    layer{end}.dim = [layer{end}.delta_order+1 1]*layer{length(layer)+layer{end}.prev}.dim(1);
end

mean_idx = length(layer);

% get the var
tmpTopology.inputDim = layer{shared_lstm_idx}.dim(1);
tmpTopology.hiddenLayerSizeLSTM = para.hiddenLayerSizeVar;
tmpTopology.usePastState = zeros(1,length(para.hiddenLayerSizeVar)); % do not use peeping hole
tmpTopology.hiddenLayerSizeFF = [];
tmpTopology.outputDim = nFreqBin*3;
tmpTopology.costFn = 'mse';
tmpTopology.LastActivation4MSE = 'linear';
layerVar = genNetworkLSTM(tmpTopology);
layerVar = layerVar(2:end-2);
layerVar{1}.prev = shared_lstm_idx - length(layer)-1;
if stage==1
    for i=1:length(layerVar)
        layerVar{i}.update = 0;
        layerVar{i}.skipBP = 1;
    end
    % always output unit variance
    layerVar{end}.W = zeros(layerVar{end}.dim(1), layerVar{end}.dim(2));
    layerVar{end}.b = zeros(layerVar{end}.dim(1),1);
    layerVar{end}.b(nFreqBin+1:nFreqBin*2) = -log(para.MSECostWeightSDA(2)^2);
    layerVar{end}.b(nFreqBin*2+1:nFreqBin*3) = -log(para.MSECostWeightSDA(3)^2);
end
layer = [layer layerVar];

layer{end+1}.name = 'exp';
layer{end}.prev = -1;
layer{end}.dim = [1 1]* layer{end-1}.dim(1);

var_idx = length(layer);

if para.useWav
    log_idx = ReturnLayerIdxByName(layer, 'log');
    if para.useCMN
        CleanLayer = layer(1:log_idx(1)+1);
    else
        CleanLayer = layer(1:log_idx(1));
    end
    CleanLayer{1}.inputIdx = 2;
    CleanLayer{1}.dim(:) = 1;
    CleanLayer{2}.dim(:) = CleanLayer{2}.dim(:)/para.nCh;
    CleanLayer{3}.dim(:) = CleanLayer{3}.dim(:)/para.nCh;
    CleanLayer(4) = [];
else
    CleanLayer{1}.name = 'Input';
    CleanLayer{end}.inputIdx = 2;
    CleanLayer{end}.dim = [1 1]*(para.fft_len/2+1);
end

CleanLayer{end+1}.name = 'Delta';
CleanLayer{end}.prev = -1;
CleanLayer{end}.delta_order = 2;
CleanLayer{end}.dim = [CleanLayer{end}.delta_order+1 1]*CleanLayer{length(CleanLayer)+CleanLayer{end}.prev}.dim(1);

layer = [layer CleanLayer];

layer{end+1}.name = 'LL_Gaussian';
layer{end}.prev = [mean_idx-length(layer) var_idx-length(layer) -1];
layer{end}.dim = [1 layer{end-1}.dim(1)];

layer = FinishLayer(layer);
end

author: singaxiong

repo: SignalGraph


text:
%DEMO_AUDIODENOISE  Audio denoising using thresholding
%
%   This demos shows how to do audio denoising using thresholding
%   of WMDCT transform.
%
%   The signal is transformed using an orthonormal WMDCT transform
%   followed by a thresholding. Then the signal is reconstructed
%   and compared with the original.
%
%   .. figure::
%
%      Denoising
%
%      The figure shows the original signal, the noisy signal and denoised
%      signals using hard and soft threshholding applied to the WMDCT of the
%      noise signal.
%

% Load audio signal
% Use the 'glockenspiel' signal.
sig=gspi;

SigLength = 2^16;
sig = sig(1:SigLength);

% Initializations
NbFreqBands = 1024;
sigma = 0.5;
Relative_Threshold = 0.1;
tau = Relative_Threshold*sigma;

% Generate window
gamma = wilorth(NbFreqBands,SigLength);

% add noise to the signal
sigma = sigma * std(sig);
nsig = sig + sigma * randn(size(sig));

% Compute wmdct coefficients
c = wmdct(nsig,gamma,NbFreqBands);

% Hard Thresholding
chard=thresh(c,tau);

% Reconstruct
hrec = real(iwmdct(chard,gamma));

% Soft thresholding
csoft=thresh(c,tau,'soft');

% Reconstruct
srec = real(iwmdct(csoft,gamma));

% Plot
figure(1);
subplot(4,1,1); plot(sig); legend('Original');
subplot(4,1,2); plot(nsig); legend('Noisy');
subplot(4,1,3); plot(hrec); legend('Hard threshold');
subplot(4,1,4); plot(srec); legend('Soft threshold');

% Results
InputSNR = 20 *log10(std(sig)/std(nsig-sig));
OutputSNR_h = 20 *log10(std(sig)/std(hrec-sig));
OutputSNR_s = 20 *log10(std(sig)/std(srec-sig));

fprintf(' RESULTS:\n');
fprintf('      Input SNR: %f dB.\n',InputSNR);
fprintf('      Output SNR (hard): %f dB.\n',OutputSNR_h);
fprintf('      Output SNR (soft): %f dB.\n',OutputSNR_s);
fprintf(' Signals are stored in variables sig, nsig, hrec, srec\n');

author: ltfat

repo: ltfat


text:
function [priormeans,posteriormeans,covmatrix,rbfposteriormeans,rbfcovmatrix] = gpnddisimPredict(model,predtimes,predict_rna,with_obsnoise);

% GPASIMPREDICT Compute predictions (means and a covariance matrix)
% of POL2 and RNA values for the GPASIM model.
%
% FORMAT
%---------------------------------
% DESC computes predictions for the asynchronous Gaussian
% process single input motif model.
%
% ARG model : the model for which the gradient is computed.
%
% ARG pol2times : the time points where predictions for POL2 are needed
%
% ARG rnatime : the time points where predictions for RNA are needed
%
% RETURN means : the predicted mean values, first the POL2
% predictions and then the RNA predictions.
%
% RETURN covmatrix : the covariance matrix between the
% predictions; for example, the diagonal values are the variances
% of each prediction.
%---------------------------------
%
% SEEALSO : gpasimCreate
%
% COPYRIGHT : Jaakko Peltonen, 2011

% GPASIMPREDICT

%predtimes
%pause

if nargin < 4,
  with_obsnoise = 1;
end

numGenes=model.numGenes;
if predict_rna==0,
  numGenes=0;
end;

% compute prior means
if iscell(predtimes)==0,
  pol2priormeans=ones(size(predtimes,1),1)*model.simMean;
else
  pol2priormeans=ones(size(predtimes{1},1),1)*model.simMean;
end;
%model.pol2mean;

if numGenes>0,
  % Mean for the mRNA is nonconstant over time and depends on the
  % B,D,S parameters and on the POL2 mean
  Bj=model.B(1);
  Dj=model.D(1);
  Sj=model.S(1);
  if model.use_disimstartmean==1,
    disimStartMean=model.disimStartMean(1);
  end;
end;

if numGenes>0,
  % compute the RNA mean curve
%  if model.use_disimstartmean==1,
%    rnapriormeans=(Bj+model.simMean*Sj)/Dj+(disimStartMean-(Bj+model.simMean*Sj)/Dj)*exp(Dj*(-predtimes));
%  else
%    rnapriormeans=(Bj+model.simMean*Sj)/Dj*ones(size(predtimes));
%  end;

  rnapriormeans=[];
  tempind1=1;
  for k=1:numGenes,
    if iscell(predtimes)==0,
      nt=length(predtimes);
    else
      nt=length(predtimes{k+1});
    end;

    rnapriormeans=[rnapriormeans;nan*ones(nt,1)];
    if (model.use_disimstartmean==1),
      if iscell(predtimes)==0,
        tempt=predtimes;
      else
        tempt=predtimes{k+1};
      end;
      delayedt=tempt-model.delay(k);

      I=find(delayedt<0);
      delayedt(I)=0;
      rnapriormeans(tempind1:tempind1+nt-1)=...
          model.disimStartMean(k)*exp(model.D(k)*(-tempt)) ...
          +(model.B(k)/model.D(k))*(1-exp(-model.D(k)*tempt)) ...
          +(model.simMean*model.S(k)/model.D(k))*(1-exp(-model.D(k)*delayedt));
    else
      if iscell(predtimes)==0,
        tempt=predtimes;
      else
        tempt=predtimes{k+1};
      end;
      delayedt=tempt-model.delay(k);
      I=find(delayedt<0);
      delayedt(I)=0;
      rnapriormeans(tempind1:tempind1+nt-1)=...
          ((model.B(k)+model.simMean*model.S(k))/model.D(k))*exp(model.D(k)*(-tempt))...
          +((model.B(k)+model.simMean*model.S(k))/model.D(k))*(1-exp(-model.D(k)*delayedt));
    end;
    tempind1=tempind1+nt;
  end;
  %size(rnapriormeans)
  %size(pol2priormeans)
end;

if 1,
if with_obsnoise,
  % This version of K_new does include observation noise
  K_new=kernCompute(model.kern, predtimes);
else
  % This version of K_new does not include observation noise
  K_new=kernCompute(model.kern, predtimes, predtimes);
end

predmodeltimes=model.t;
if (iscell(predtimes)==1) && (iscell(model.t)==0),
  predmodeltimes={model.t,model.t};
end;
if (iscell(predtimes)==0) && (iscell(model.t)==1),
  predtimes={predtimes,predtimes};
end;

K_new_old=kernCompute(model.kern, predtimes, predmodeltimes);
K_old=model.K;
K_old_new=K_new_old';
end;

%K_old
%pause

if 0,
K_old_ndsim=ndsimKernCompute(model.kern.comp{1}.comp{1},model.t);
K_old_nddisim=nddisimKernCompute(model.kern.comp{1}.comp{2},model.t);
K_old_nddisimXndsim=nddisimXndsimKernCompute(model.kern.comp{1}.comp{2},model.kern.comp{1}.comp{1},model.t);
K_old=[K_old_ndsim K_old_nddisimXndsim';K_old_nddisimXndsim K_old_nddisim];
K_old=real(K_old);

K_new_ndsim=ndsimKernCompute(model.kern.comp{1}.comp{1},predtimes);
K_new_nddisim=nddisimKernCompute(model.kern.comp{1}.comp{2},predtimes);
K_new_nddisimXndsim=nddisimXndsimKernCompute(model.kern.comp{1}.comp{2},model.kern.comp{1}.comp{1},predtimes);
K_new=[K_new_ndsim K_new_nddisimXndsim';K_new_nddisimXndsim K_new_nddisim];
K_new=real(K_new);

K_new_old_ndsim=ndsimKernCompute(model.kern.comp{1}.comp{1},predtimes,model.t);
K_new_old_nddisim=nddisimKernCompute(model.kern.comp{1}.comp{2},predtimes,model.t);
K_new_old_ndsimXnddisim=nddisimXndsimKernCompute(model.kern.comp{1}.comp{2},model.kern.comp{1}.comp{1},model.t,predtimes)';
K_new_old_nddisimXndsim=nddisimXndsimKernCompute(model.kern.comp{1}.comp{2},model.kern.comp{1}.comp{1},predtimes,model.t);
K_new_old=[K_new_old_ndsim K_new_old_ndsimXnddisim;K_new_old_nddisimXndsim K_new_old_nddisim];
K_new_old=real(K_new_old);

K_old_new=K_new_old';

noisekern_new=kernCompute(model.kern.comp{2}, predtimes);
noisekern_old=kernCompute(model.kern.comp{2}, model.t);
K_old=K_old+noisekern_old;
K_new=K_new+noisekern_new;
end;

tempm=model.m;

% If we do not want to predict using RNA observations, throw out kernel parts related to RNA
if numGenes==0,
  if iscell(model.t)==0,
    ot=length(model.t);
    nt=length(predtimes);
  else
    % ot=0;
    % nt=0;
    % for k=1:length(model.t),
    %   ot=ot+size(model.t{k},1);
    %   nt=nt+size(predtimes{k},1);
    % end;
    ot=length(model.t{1});
    nt=length(predtimes{1});
  end;
  K_old=K_old(1:ot,1:ot);
  K_new_old=K_new_old(1:nt,1:ot);
  K_old_new=K_old_new(1:ot,1:nt);
  K_new=K_new(1:nt,1:nt);
  tempm=tempm(1:ot);
end;

if numGenes>0,
  priormeans=[pol2priormeans;rnapriormeans];
else
  priormeans=pol2priormeans;
end;

%predict_rna
%K_old
%K_new_old
%pause
%figure; imagesc(K_old);
%pause

posteriormeans=priormeans+K_new_old*(K_old\tempm);
%covmatrix=K_new-K_new_old*inv(K_old)*K_old_new;

covmatrix=K_new-K_new_old*(K_old\K_old_new);
if(min(diag(covmatrix))<0),
  % Try omitting the first row and column of the kernel
  % since it might be all zeroes

  if iscell(model.t)==0,
    ot_pol2=length(model.t);
    ot_rna=length(model.t);
    nt_pol2=length(predtimes);
    nt_rna=length(predtimes);
  else
    ot_pol2=size(model.t{1},1);
    ot_rna=size(model.t{2},1);
    nt_pol2=size(predtimes{1},1);
    nt_rna=size(predtimes{2},1);
  end;

  if (numGenes>0),
    okentries=[(2:ot_pol2) (ot_pol2+2:ot_pol2+ot_rna)];
  else
    okentries=[2:ot_pol2];
  end;

  K_old=K_old(okentries,okentries);
  K_new_old=K_new_old(:,okentries);
  K_old_new=K_old_new(okentries,:);
  covmatrix=K_new-K_new_old*(K_old\K_old_new);
end;

if 0,
  % Compute predictions for the driving RBF-kernel GP; only valid for
  % the SIM-DISIM model.

  % extract RBF transforms settings from DISIM part of the SIM-DISIM kernel
  disimkern=model.kern.comp{1}.comp{2};
  if isfield(disimkern, 'options') ...
        && isfield(disimkern.options, 'isNegativeS') ...
        && kern.options.isNegativeS,
    error('DISIM kern uses negative variance, cannot decide transform settings for RBF-SIM-DISIM prediction');
  else
    rbftransformsettings{1}=disimkern.transforms(1).transformsettings; % setting for inverse width
    rbftransformsettings{2}=[0 10]; % setting for RBF variance, not
                                    % used in GPNDDISIM model
  end

  % create a RBF kernel to match the settings in the NDSIM-NDDISIM model
  rbfkern=kernCreate(predtimes,'rbf');
  rbfkern=kernExpandParamTransformSettings(rbfkern,rbftransformsettings);
  %[rbfpars,rbfnams]=kernExtractParam(rbfkern);
  [modelpars,modelnams]=gpnddisimExtractParam(model);
 % modelpars
 % modelnams
 % pause
  fprintf(1,'Calling gpdisimExpandParam\n');
  model=gpnddisimExpandParam(model,modelpars);
  fprintf(1,'Calling gpdisimExpandParam done\n');

  rbfpar_inversewidth=modelpars(1);  % assumes that inverse width is parameter 1 in the NDSIM-NDDISIM model
  rbfpar_variance=sigmoidabTransform(1, 'xtoa', rbftransformsettings{2});
  rbfpars(1)=rbfpar_inversewidth;
  rbfpars(2)=rbfpar_variance;
  rbfkern=kernExpandParam(rbfkern,rbfpars);
%  rbfkern
%  rbftransformsettings

  % create fake SIM and DISIM kernels corresponding to the NDSIM
  % and NDDISIM kernels... this is just a quick hack

  K_rbf=kernCompute(rbfkern,predtimes);
  K_rbf_sim=simXrbfKernCompute(model.kern.comp{1}.comp{1},rbfkern,model.t,predtimes)';
  K_rbf_disim=disimXrbfKernCompute(model.kern.comp{1}.comp{2},rbfkern,model.t,predtimes)';
  K_rbf_old=[K_rbf_sim K_rbf_disim];
  K_old_rbf=K_rbf_old';
  %size(K_rbf_old)
  %size(K_old)
  %size(model.m)

  K_old=model.K;
  rbfposteriormeans=0+K_rbf_old*(K_old\model.m);
  rbfcovmatrix=K_rbf-K_rbf_old*inv(K_old)*K_old_rbf;

%  rbfposteriormeans=real(rbfposteriormeans);
%  rbfcovmatrix=real(rbfcovmatrix);
else
  rbfposteriormeans=[];
  rbfcovmatrix=[];
end;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function area = hypersphere_01_area ( dim_num )

%*****************************************************************************80
%
%% HYPERSPHERE_01_AREA computes the surface area of a unit sphere in ND.
%
%  Discussion:
%
%    The unit sphere in ND satisfies:
%
%      sum ( 1 <= I <= DIM_NUM ) X(I) * X(I) = 1
%
%    Results for the first few values of N are:
%
%    DIM_NUM   Area
%
%     2    2        * PI
%     3    4        * PI
%     4  ( 2 /   1) * PI^2
%     5  ( 8 /   3) * PI^2
%     6  ( 1 /   1) * PI^3
%     7  (16 /  15) * PI^3
%     8  ( 1 /   3) * PI^4
%     9  (32 / 105) * PI^4
%    10  ( 1 /  12) * PI^5
%
%    For the unit sphere, Area(DIM_NUM) = DIM_NUM * Volume(DIM_NUM)
%
%    Sphere_Unit_Area ( DIM_NUM ) = 2 * PI^(DIM_NUM/2) / Gamma ( DIM_NUM / 2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the space.
%
%    Output, real SPHERE_UNIT_AREA_ND, the area of the sphere.
%
  if ( mod ( dim_num, 2 ) == 0 )
    m = floor ( dim_num / 2 );
    area = 2.0 * pi^m;
    for i = 1 : m-1
      area = area / i;
    end
  else
    m = floor ( ( dim_num - 1 ) / 2 );
    area = pi^m * 2.0^dim_num;
    for i = m+1 : 2*m
      area = area / i;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [cs,mineral] = loadCIF(fname,varargin)
% import crystal symmetry from cif file
%
% if cif file not found and input name is a valid COD entry, this function
% tries to download the file from <http://www.crystallography.net/cif/
% http://www.crystallography.net/cif/>
%
% Syntax
%   loadCIF('5000035.cif')
%   loadCIF(5000035)       % lookup online
%
% See also
% symmetry

if ~iscell(fname)
  if isnumeric(fname)
    fname = copyonline(fname);
  end
  [pathstr, name, ext] = fileparts(char(fname));

  if isempty(ext), ext = '.cif';end
  if isempty(pathstr) && ~exist([name,ext],'file')
    pathstr = mtexCifPath;
  end

  % load file
  if ~exist(fullfile(pathstr,[name ext]),'file')
    try
      fname = copyonline(fname);
    catch %#ok<CTCH>
      dir(fullfile(mtexCifPath,'*.cif'))
      error('I could not find the corresponding cif. Above you see the list of localy avaible cif files.')
    end
  else
    fname = fullfile(pathstr,[name ext]);
  end
  str = file2cell(fname);
else
  str = fname;
  name = '';
end

try
  % get a name for it
  mineral_names = {...
    '_chemical_name_mineral',...
    '_chemical_name_systematic',...
    '_chemical_formula_structural',...
    '_chemical_formula_sum'};

  for alias = mineral_names
    mineral = extract_token(str,alias{:});
    if ~isempty(mineral), break; end
  end

  % find space group
  group_aliases = {...
    '_symmetry_space_group_name_H-M',...
    '_symmetry_point_group_name_H-M',...
    '_symmetry_cell_setting'};
  for gp = group_aliases
    group = extract_token(str,gp{:});
    if ~isempty(group), break; end
  end

  % find a,b,c
  axis = [extract_token(str,'_cell_length_a',true) ...
    extract_token(str,'_cell_length_b',true) ...
    extract_token(str,'_cell_length_c',true)];

  % find alpha, beta, gamma
  angles = [extract_token(str,'_cell_angle_alpha',true) ...
    extract_token(str,'_cell_angle_beta',true) ...
    extract_token(str,'_cell_angle_gamma',true)];

  if length(axis)<3,
%     warning('crystallographic axis mismatch');
    axis = [1 1 1];
  end
  if length(angles)<3,
%     warning('crystallographic angles mismatch');
    angles = [90 90 90];
  end

  assert(~isempty(group));

  cs = crystalSymmetry(group,axis,angles*degree,'mineral',mineral);

catch
  error(['Error reading cif file', fname]);
end

function t = extract_token(str,token,numeric)

pos = strmatch(token,str);
if ~isempty(pos)
  t = strtrim(regexprep(str{pos(1)},[token '|'''],''));
  if ~isempty(t)
    if nargin>2 && numeric
      t = sscanf(t,'%f');
    end
  end
else
  t = '';
end

function fname = copyonline(cod)

try
  if isnumeric(cod)
    cod = num2str(cod);
  end

  fname = fullfile(mtexCifPath,[cod '.cif']);
  if exist(fname,'file')
    return
  end

  if ~isempty(str2num(cod))
    disp('CIF-File from Crystallography Open Database')
    disp(['> download : http://www.crystallography.net/cif/' cod '.cif'])
    cif = urlread(['http://www.crystallography.net/cif/' cod '.cif']);
  else
    cif = urlread(cod);
  end
catch
  disp('> unluckily failed to find cif-file')
  error('CIF-file not valid online')
end

fname = fullfile(mtexCifPath,[cod '.cif']);

fid = fopen(fname,'w');
fwrite(fid,cif);
fclose(fid);
disp(['> copied to: ' fname]);

author: mtex-toolbox

repo: mtex

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


text:
% sbplot() - create axes in arbitrary subplot grid positions and sizes
%
% Usage:  >> axis_handle = sbplot(v,h,index)
%         >> axis_handle = sbplot(v,h,[index1 index2])
%         >> axis_handle = sbplot(v,h,[index1 index2],axprop,..)
%         >> axis_handle = sbplot(v,h,[index1 index2],'ax',handle,axprop,..)
%
% Inputs:
%   v,h    - Integers giving the vertical and horizontal ranks of the tiling.
%   index  - Either a single subplot index, in which case the command
%            is equivalent to subplot, or a two-element vector giving
%            the indices of two corners of the sbplot() area according
%            to subplot() convention (e.g., left-to-right, top-to-bottom).
%   axprop - Any axes property(s), e.g., >> sbplot(3,3,3,'color','w')
%   handle - Following keyword 'ax', sbplot tiles the given axes handle
%            instead of the whole figure
%
% Output:
%  axis_handle - matlab axis handle
%
% Note:
%   sbplot is essentially the same as the subplot command except that
%   sbplot axes may span multiple tiles. Also, sbplot() will not erase
%   underlying axes.
%
%  Examples: >> sbplot(3,3,6);plot(rand(1,10),'g');
%            >> sbplot(3,3,[7 2]);plot(rand(1,10),'r');
%            >> sbplot(8,7,47);plot(rand(1,10),'b');
%
% Authors: Colin Humphries, Arnaud Delorme & Scott Makeig, SCCN/INC/UCSD, La Jolla, June, 1998

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

% reformatted by Scott Makeig, 6/10/98
% 12/22/00 test nargin<3 -sm
% 01/21/01 added (recursive) axes option 'ax' -sm
% 01-25-02 reformated help & licence -ad

function [out] = sbplot(m,n,gridpos,varargin) % varargin is std. matlab arg list

if nargin<3
   error('            requires >=3 arguments');
end

if nargin>3 & strcmp(varargin{1},'ax')
   pos = get(varargin{2},'Position'); %  sbplot(3,1,[2 3]) -> 0.4111 0.1100 0.4939 0.815
   varargin = {varargin{3:end}};
else
   pos = get(gcf,'DefaultAxesPosition'); %  [0.1300    0.1100    0.7750    0.815]
end

Xpad = pos(1);      % lower-left distance from left side of figure
Ypad = pos(2);      % lower-left distance from bottom of figure

Xlen = pos(3);      % axes width
Ylen = pos(4);      % axes height

if n == 2
  xspace = Xlen*0.27/(n-0.27);      % xspace between axes as per subplot
else
  xspace = (0.9*Xlen)*0.27/(n-0.9*0.27);
end

if m == 2
  yspace = Ylen*0.27/(m-0.27);      % yspace between axes as per subplot
else                                % WHY Xlen (.775) instead of Ylen (.815) ??
  yspace = (0.9*Ylen)*0.27/(m-0.9*0.27);
end

xlength = (Xlen-xspace*(n-1))/n;  % axes width
ylength = (Ylen-yspace*(m-1))/m;  % axes height

% Convert tile indices to grid positions
if length(gridpos) == 1
  xgridpos(1) = mod(gridpos,n);   % grid position
  if xgridpos(1) == 0
    xgridpos(1) = n;
  end
  xgridpos(2) = 1;                   % grid length
  ygridpos(1) = m-ceil(gridpos/n)+1; % grid position
  ygridpos(2) = 1;                   % grid length
else
  xgridpos(1) = mod(gridpos(1),n);
  if xgridpos(1) == 0
    xgridpos(1) = n;
  end
  tmp = mod(gridpos(2),n);
  if tmp == 0
    tmp = n;
  end
  if tmp > xgridpos(1)
    xgridpos(2) = tmp-xgridpos(1)+1;
  else
    xgridpos(2) = xgridpos(1)-tmp+1;
    xgridpos(1) = tmp;
  end

  ygridpos(1) = m-ceil(gridpos(1)/n)+1;
  tmp = m-ceil(gridpos(2)/n)+1;
  if tmp > ygridpos(1)
    ygridpos(2) = tmp-ygridpos(1)+1;
  else
    ygridpos(2) = ygridpos(1)-tmp+1;
    ygridpos(1) = tmp;
  end
end

% Calculate axes coordinates
position(1) = Xpad+xspace*(xgridpos(1)-1)+xlength*(xgridpos(1)-1);
position(2) = Ypad+yspace*(ygridpos(1)-1)+ylength*(ygridpos(1)-1)-0.03;
position(3) = xspace*(xgridpos(2)-1)+xlength*xgridpos(2);
position(4) = yspace*(ygridpos(2)-1)+ylength*ygridpos(2);

% Create new axes
ax = axes('Position',position,varargin{:});

% Output axes handle
if nargout > 0
  out = ax;
end

author: PatternRecognition

repo: OpenBMI


text:
% Prim's minimal spanning tree algorithm
% Prim's alg idea:
%  start at any node, find closest neighbor and mark edges
%  for all remaining nodes, find closest to previous cluster, mark edge
%  continue until no nodes remain
%
% INPUTS: graph defined by adjacency matrix, nxn
% OUTPUTS: matrix specifying minimum spanning tree (subgraph), nxn
%
% Other routines used: isConnected.m
% GB: Oct 7, 2012

function tr = minSpanTree(adj)

% check if graph is connected:
if not(isConnected(adj)); printf('This graph is not connected. No spanning tree exists.\n'); return; end

n = length(adj); % number of nodes
tr = zeros(n);   % initialize tree

adj(find(adj==0))=inf; % set all zeros in the matrix to inf

conn_nodes = 1;        % nodes part of the min-span-tree
rem_nodes = [2:n];     % remaining nodes

while length(rem_nodes)>0

  [minlink]=min(min(adj(conn_nodes,rem_nodes)));
  ind=find(adj(conn_nodes,rem_nodes)==minlink);

  [ind_i,ind_j] = ind2sub([length(conn_nodes),length(rem_nodes)],ind(1));

  i=conn_nodes(ind_i); j=rem_nodes(ind_j); % gets back to adj indices
  tr(i,j)=1; tr(j,i)=1;
  conn_nodes = [conn_nodes j];
  rem_nodes = setdiff(rem_nodes,j);

end

author: aeolianine

repo: octave-networks-toolbox


text:
function e = polyterm_exponent ( action, e )

%*****************************************************************************80
%
%% POLYTERM_EXPONENT gets or sets the exponents for the polynomial term.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION.
%    'GET' asks the routine to return the current values in E.
%    'SET' asks the routine to set the current values to E.
%
%    Input/output, integer E(3), storage used to set or get values.
%
  persistent e_save

  if ( s_eqi ( action, 'GET' ) )

    e(1:3) = e_save(1:3);

  elseif ( s_eqi ( action, 'PRINT' ) )

    fprintf ( 1, '\n' );

    if ( all ( e_save(1:3) == 0 ) )

      fprintf ( 1, 'P(X,Y,Z) = 1' );

    else

      fprintf ( 'P(X,Y,Z) = ' );

      if ( e_save(1) == 0 )

      elseif ( e_save(1) == 1 )

        fprintf ( 1, 'X ' );

      else

        fprintf ( 1, 'X^%d', e_save(1) );

      end

      if ( e_save(2) == 0 )

      elseif ( e_save(2) == 1 )

        fprintf ( 1, ' Y ' );

      else

        fprintf ( 1, ' Y^%d', e_save(2) );

      end

      if ( e_save(3) == 0 )

      elseif ( e_save(3) == 1 )

        fprintf ( 1, ' Z' );

      else

        fprintf ( 1, ' Z^%d', e_save(3) );

      end

    end

    fprintf ( 1, '\n' );

  elseif ( s_eqi ( action, 'SET' ) )

    e_save(1:3) = e(1:3);

  end

  return
end


text:
function Y = scale(X,S,dims)
%SCALE Scale along specified dimensions of tensor.
%
%   Y = SCALE(X,S,DIMS) scales the tensor X along the dimension(s)
%   specified in DIMS using the scaling data in S. If DIMS contains
%   only one dimension, then S can be a column vector. Otherwise, S
%   should be a tensor.
%
%   Examples
%   X = tenones([3,4,5]);
%   S = 10 * [1:5]'; Y = scale(X,S,3)
%   S = tensor(10 * [1:5]',5); Y = scale(X,S,3)
%   S = tensor(1:12,[3 4]); Y = scale(X,S,[1 2])
%   S = tensor(1:12,[3 4]); Y = scale(X,S,-3)
%   S = tensor(1:60,[3 4 5]); Y = scale(X,S,1:3)
%
%   See also TENSOR, TENSOR/COLLAPSE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

dims = tt_dimscheck(dims,ndims(X));
remdims = setdiff(1:ndims(X),dims);

% Convert to a matrix so that each column of A can be scaled by a
% vectorized version of S.
A = double(tenmat(X,dims,remdims));

switch(class(S))
    case {'tensor'}
        if ~isequal(size(S), X.size(dims))
            error 'Size mismatch';
        end
        % Vectorize S.
        S = double(tenmat(S,1:ndims(S),[]));
    case {'double'}
        if size(S,1) ~= X.size(dims)
            error 'Size mismatch';
        end
    otherwise
        error('Invalid scaling factor');
end

[m,n] = size(A);

% If the size of S is pretty small, we can convert it to a diagonal matrix
% and multiply by A. Otherwise, we scale A column-by-column.
if (m <= n)
    B = diag(S) * A;
else
    B = zeros(size(A));
    for j = 1:n
        B(:,j) = S .* A(:,j);
    end
end

% Convert the matrix B back into a tensor and return.
Y = tensor(tenmat(B,dims,remdims,X.size));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm