text:
%PoseGraph Pose graph

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
classdef PoseGraph < handle

    properties
        graph

        ngrid
        center
        cellsize
    end

    methods

        function pg = PoseGraph(filename, varargin)
            % parse the file data
            % we assume g2o format
            %    VERTEX* vertex_id X Y THETA
            %    EDGE* startvertex_id endvertex_id X Y THETA IXX IXY IYY IXT IYT ITT
            % vertex numbers start at 0

            opt.laser = false;

            opt = tb_optparse(opt, varargin);

            pg.graph = PGraph(3, 'distance', 'SE2');

            fp = fopen(filename, 'r');
            assert(fp > 0, 'Can''t open file %s', filename);

            toroformat = false;
            nlaser = 0;

            % indices into ROBOTLASER1 record for the 3x3 info matrix in column major
            % order
            g2o =  [6  7  8  7  9  10   8  10  11];
            toro = [6  7 10  7  8  11  10  11   9];

            % we keep an array pgi = vindex(gi) to map g2o vertex index to PGraph vertex index

            tic
            while ~feof(fp)

                line = fgets(fp);

                % is it a comment?
                if line(1) == '#'
                    continue;
                end

                % get keyword
                k = strfind(line, ' ');

                % and deal with it
                switch line(1:k-1)
                    case 'VERTEX_SE2'
                        % g2o format vertex
                        vertex = sscanf(line(k+1:end), '%d %f %f %f')';
                        v = pg.graph.add_node(vertex(2:4));
                        vindex(vertex(1)+1) = v;
                        vd.type = 'vertex';
                        pg.graph.setvdata(v, vd);

                    case 'VERTEX_XY'
                        vertex = sscanf(line(k+1:end), '%d %f %f')';
                        v = pg.graph.add_node(vertex(2:4));
                        vindex(vertex(1)+1) = v;
                        vd.type = 'landmark';
                        pg.graph.setvdata(v, vd);

                    case 'EDGE_SE2'
                        % g2o format edge
                        edge = sscanf(line(k+1:end), '%f')';
                        v1 = vindex(edge(1)+1);
                        v2 = vindex(edge(2)+1);

                        % create the edge
                        e = pg.graph.add_edge(v1, v2);

                        % create the edge data as a structure
                        %  X  Y  T
                        %  3  4  5
                        ed.mean = edge(3:5);

                        % IXX IXY IXT IYY IYT ITT
                        %   6   7   8   9  10  11
                        ed.info = reshape(edge(g2o), [3 3]);

                        % and attach it
                        pg.graph.setedata(e, ed);

                    case 'VERTEX2'
                        toroformat = true;
                        vertex = sscanf(line(k+1:end), '%d %f %f %f')';
                        v = pg.graph.add_node(vertex(2:4));
                        vindex(vertex(1)+1) = v;
                        vd.type = 'vertex';
                        pg.graph.setvdata(v, vd);

                    case 'EDGE2'
                        toroformat = true;
                        edge = sscanf(line(k+1:end), '%f')';
                        v1 = vindex(edge(1)+1);
                        v2 = vindex(edge(2)+1);

                        % create the edge
                        e = pg.graph.add_edge(v1, v2);
                        % create the edge data as a structure
                        %  X  Y  T
                        %  3  4  5
                        ed.mean = edge(3:5);

                        % IXX IXY IYY ITT IXT IYT
                        %   6   7   8   9  10  11
                        ed.info = reshape(edge(toro), [3 3]);
                        % and attach it
                        pg.graph.setedata(e, ed);

                    case 'ROBOTLASER1'
                        if ~opt.laser
                            continue;
                        end

                        % laser records are associated with the immediately preceding VERTEX record
                        [laser,n] = sscanf(line(k+1:end), '%f');
                        nbeams = laser(8);
                        vd.theta = [0:nbeams-1] * laser(4) + laser(2);
                        vd.range = laser(9:8+nbeams)';
                        vd.time = laser(21+nbeams);
                        pg.graph.setvdata(v, vd);
                        nlaser = nlaser + 1;

                    otherwise
                        error('RTB:posegraph:badfile', 'Unexpected line  <%s> in %s', line(1:k-1), filename);
                end
            end
            elapsed = toc;

            fclose(fp);

            if toroformat
                fprintf('loaded TORO/LAGO format file: %d nodes, %d edges in %.2f sec\n', pg.graph.n, pg.graph.ne, elapsed);
            else
                fprintf('loaded g2o format file: %d nodes, %d edges in %.2f sec\n', pg.graph.n, pg.graph.ne, elapsed);
                if nlaser > 0
                    fprintf('  %d laser scans: %d beams, fov %g to %g deg, max range %g\n', ...
                        nlaser, nbeams, [laser(2) sum(laser(2:3))]*180/pi, laser(5) );
                end
            end
        end

        function [r, theta] = scan(pg, n)
            vd = pg.graph.vdata(n);
            r = vd.range;
            theta = vd.theta;
        end

        function [X,Y] = scanxy(pg, n)
            vd = pg.graph.vdata(n);

            [x,y] = pol2cart(vd.theta, vd.range);
            if nargout == 1
                X = [x; y];
            elseif nargout == 2
                X = x; Y = y;
            end
        end

        function plot_scan(pg, n)
            for i=n(:)'
                [x,y] = pg.scanxy(i);
                plot(x, y, '.', 'MarkerSize', 10);
                pause
            end
        end

        function xyt = pose(pg, i)
            xyt = pg.graph.coord(i);
        end

        function t = time(pg, n)
            t = pg.graph.vdata(n).time;
        end

        function plot(pg, varargin)
            pg.graph.plot(varargin{:});
            xlabel('x')
            ylabel('y')
            grid on
        end

        function world = scanmap(pg, varargin)

            opt.center = [75 50];
            opt.ngrid = 3000;
            opt.cellsize = 0.1;

            pg = tb_optparse(opt, varargin, pg);

            h = waitbar(0, 'rendering a map');

            world = zeros(pg.ngrid, pg.ngrid, 'int32');
            for i=1:1:pg.graph.n

                if rem(i, 20) == 0
                    waitbar(i/pg.graph.n, h)
                end

                xy = pg.scanxy(i);
                [r,theta] = pg.scan(i);
                xy(:,r>40) = [];
                xyt = pg.graph.coord(i);

                xy = SE2(xyt) * xy;

                % start of each ray
                [x1,y1] = pg.w2g(xyt(1:2));

                for s=1:numcols(xy)

                    % end of each ray
                    [x2,y2] = pg.w2g(xy(:,s));

                    % all cells along the ray
                    p = bresenham(x1, y1, x2, y2);
                    try
                        k = sub2ind(size(world), p(:,1), p(:,2));

                        k1 = k(1:end-1); k2 = k(end);
                        world(k1) = world(k1) - 1;
                        world(k2) = world(k2) + 1;
                    catch me
                        % come here if any point on the ray is outside the grid
                        % silently ignore it
                    end
                end
            end
            close(h)
            %idisp(world)
        end

        function [gx,gy] = w2g(pg, w)
            dd = 0.10;

            w = w(:) + pg.center(:);
            g = round(w/pg.cellsize);
            gx = g(1); gy = g(2);
        end

        function plot_occgrid(pg, w)

            x = [1:numcols(w)]*pg.cellsize - pg.center(1);
            y = [1:numrows(w)]*pg.cellsize - pg.center(2);
            w(w<0) = -1;
            w(w>0) = 1;
            w=-w;
            idisp(w, 'nogui', 'xydata', {x, y})
            xlabel('x'); ylabel('y');
        end

        %   This source code is part of the graph optimization package
        %   deveoped for the lectures of robotics2 at the University of Freiburg.
        %
        %     Copyright (c) 2007 Giorgio Grisetti, Gian Diego Tipaldi
        %
        %   It is licences under the Common Creative License,
        %   Attribution-NonCommercial-ShareAlike 3.0
        %
        %   You are free:
        %     - to Share - to copy, distribute and transmit the work
        %     - to Remix - to adapt the work
        %
        %   Under the following conditions:
        %
        %     - Attribution. You must attribute the work in the manner specified
        %       by the author or licensor (but not in any way that suggests that
        %       they endorse you or your use of the work).
        %
        %     - Noncommercial. You may not use this work for commercial purposes.
        %
        %     - Share Alike. If you alter, transform, or build upon this work,
        %       you may distribute the resulting work only under the same or
        %       similar license to this one.
        %
        %   Any of the above conditions can be waived if you get permission
        %   from the copyright holder.  Nothing in this license impairs or
        %   restricts the author's moral rights.
        %
        %   This software 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.

        %ls-slam.m
        %this file is released under the creative common license

        %solves a graph-based slam problem via least squares
        %vmeans: matrix containing the column vectors of the poses of the vertices
        %	 the vertices are odrered such that vmeans[i] corresponds to the ith id
        %eids:	 matrix containing the column vectors [idFrom, idTo]' of the ids of the vertices
        %	 eids[k] corresponds to emeans[k] and einfs[k].
        %emeans: matrix containing the column vectors of the poses of the edges
        %einfs:  3d matrix containing the information matrices of the edges
        %	 einfs(:,:,k) refers to the information matrix of the k-th edge.
        %n:	 number of iterations
        %newmeans: matrix containing the column vectors of the updated vertices positions

        function g2 = optimize(pg, varargin)

            opt.iterations = 10;
            opt.animate = false;
            opt.retain = false;

            opt = tb_optparse(opt, varargin);

            g2 = PGraph(pg.graph);  % deep copy

            eprev = Inf;
            for i=1:opt.iterations
                if opt.animate
                    if ~opt.retain
                        clf
                    end
                    g2.plot();
                    pause(0.5)
                end

                [vmeans,energy] = linearize_and_solve(g2);
                g2.setcoord(vmeans);

                if energy >= eprev
                    break;
                end
                eprev = energy;
            end;

            pg.graph = g2;
        end
    end % methods
end % classdef

%computes the taylor expansion of the error function of the k_th edge
%vmeans: vertices positions
%eids:   edge ids
%emeans: edge means
%k:	 edge number
%e:	 e_k(x)
%A:	 d e_k(x) / d(x_i)
%B:	 d e_k(x) / d(x_j)
%function [e, A, B]=linear_factors(vmeans, eids, emeans, k)
function [e, A, B]=linear_factors(g, edge)
    %extract the ids of the vertices connected by the kth edge
    % 	id_i=eids(1,k);
    % 	id_j=eids(2,k);
    %extract the poses of the vertices and the mean of the edge
    %     v_i=vmeans(:,id_i);
    %     v_j=vmeans(:,id_j);
    %     z_ij=emeans(:,k);

    v = g.vertices(edge);
    v_i = g.coord(v(1));
    v_j = g.coord(v(2));
    z_ij = g.edata(edge).mean;

    %compute the homoeneous transforms of the previous solutions
    zt_ij=v2t(z_ij);
    vt_i=v2t(v_i);
    vt_j=v2t(v_j);

    %compute the displacement between x_i and x_j

    f_ij=(inv(vt_i)*vt_j);

    %this below is too long to explain, to understand it derive it by hand
    theta_i=v_i(3);
    ti=v_i(1:2,1);
    tj=v_j(1:2,1);
    dt_ij=tj-ti;

    si=sin(theta_i);
    ci=cos(theta_i);

    A= [-ci, -si, [-si, ci]*dt_ij; si, -ci, [-ci, -si]*dt_ij; 0, 0, -1 ];
    B =[  ci, si, 0           ; -si, ci, 0            ; 0, 0, 1 ];

    ztinv=inv(zt_ij);
    e=t2v(ztinv*f_ij);
    ztinv(1:2,3) = 0;
    A=ztinv*A;
    B=ztinv*B;

    % 	%compute the homogeneous transforms of the previous solutions
    % 	zt_ij=v2t(z_ij);
    % 	vt_i=v2t(v_i);
    % 	vt_j=v2t(v_j);
    % %     zt_ij = SE2(z_ij);
    % % 	vt_i = SE2(v_i);
    % % 	vt_j = SE2(v_j);
    %
    % 	%compute the displacement between x_i and x_j
    % 	%f_ij=(inverse(vt_i)*vt_j);
    % 	f_ij = vt_i.inv * vt_j;
    %
    % 	%this below is too long to explain, to understand it derive it by hand
    %       	theta_i=v_i(3);
    % 	ti=v_i(1:2);
    % 	tj=v_j(1:2);
    %       	dt_ij=tj-ti;
    %
    % 	si=sin(theta_i);
    % 	ci=cos(theta_i);
    %
    % 	A= [-ci, -si, [-si, ci]*dt_ij; si, -ci, [-ci, -si]*dt_ij; 0, 0, -1 ];
    % 	B =[  ci, si, 0           ; -si, ci, 0            ; 0, 0, 1 ];
    %
    % 	ztinv = inv(zt_ij);
    % 	e = xyt(ztinv*f_ij);
    % 	ztinv.t = 0;
    % 	A = ztinv*A;
    % 	B = ztinv*B;
end

%linearizes and solves one time the ls-slam problem specified by the input
%vmeans:   vertices positions at the linearization point
%eids:     edge ids
%emeans:   edge means
%einfs:    edge information matrices
%newmeans: new solution computed from the initial guess in vmeans
function [newmeans,energy] = linearize_and_solve(g)
    tic
    fprintf('solving');

    % H and b are respectively the system matrix and the system vector
    H=zeros(g.n*3,g.n*3);
    b=zeros(g.n*3,1);
    % this loop constructs the global system by accumulating in H and b the contributions
    % of all edges (see lecture)
    %for k=1:size(eids,2)
    fprintf('.');

    etotal = 0;
    for edge = 1:g.ne

        [e, A, B]=linear_factors(g, edge);
        omega = g.edata(edge).info;
        %compute the blocks of H^k

        % not quite sure whey SE3 is being transposed, what does that mean?
        b_i = -A'*omega*e;
        b_j = -B'*omega*e;
        H_ii = A'*omega*A;
        H_ij = A'*omega*B;
        H_jj = B'*omega*B;

        v = g.vertices(edge);
        id_i = v(1); id_j = v(2);
        %accumulate the blocks in H and b
        H((id_i-1)*3+1:id_i*3,(id_i-1)*3+1:id_i*3) = H((id_i-1)*3+1:id_i*3,(id_i-1)*3+1:id_i*3) + H_ii;
        H((id_j-1)*3+1:id_j*3,(id_j-1)*3+1:id_j*3) = H((id_j-1)*3+1:id_j*3,(id_j-1)*3+1:id_j*3) + H_jj;
        H((id_i-1)*3+1:id_i*3,(id_j-1)*3+1:id_j*3) = H((id_i-1)*3+1:id_i*3,(id_j-1)*3+1:id_j*3) + H_ij;
        H((id_j-1)*3+1:id_j*3,(id_i-1)*3+1:id_i*3) = H((id_j-1)*3+1:id_j*3,(id_i-1)*3+1:id_i*3) + H_ij';

        b((id_i-1)*3+1:id_i*3,1) = b((id_i-1)*3+1:id_i*3,1) + b_i;
        b((id_j-1)*3+1:id_j*3,1) = b((id_j-1)*3+1:id_j*3,1) + b_j;

        %NOTE on Matlab compatibility: note that we use the += operator which is octave specific
        %using H=H+.... results in a tremendous overhead since the matrix would be entirely copied every time
        %and the matrix is huge
        etotal = etotal + e'*e;
    end;
    fprintf('.');

    %note that the system (H b) is obtained only from
    %relative constraints. H is not full rank.
    %we solve the problem by anchoring the position of
    %the the first vertex.
    %this can be expressed by adding the equation
    %  deltax(1:3,1)=0;
    %which is equivalent to the following
    H(1:3,1:3) = H(1:3,1:3) + eye(3);

    SH=sparse(H);
    fprintf('.');
    deltax=SH\b;
    fprintf('.');

    %split the increments in nice 3x1 vectors and sum them up to the original matrix
    newmeans = g.coord()+reshape(deltax,3,g.n);

    %normalize the angles between -PI and PI
    for (i=1:size(newmeans,2))
        s=sin(newmeans(3,i));
        c=cos(newmeans(3,i));
        newmeans(3,i)=atan2(s,c);
    end
    dt = toc;
    fprintf('done in %.2g sec.  Total cost %g \n', dt, etotal);
    if nargout > 1
        energy = etotal;
    end
end

%   This source code is part of the graph optimization package
%   deveoped for the lectures of robotics2 at the University of Freiburg.
%
%     Copyright (c) 2007 Giorgio Grisetti, Gian Diego Tipaldi
%
%   It is licences under the Common Creative License,
%   Attribution-NonCommercial-ShareAlike 3.0
%
%   You are free:
%     - to Share - to copy, distribute and transmit the work
%     - to Remix - to adapt the work
%
%   Under the following conditions:
%
%     - Attribution. You must attribute the work in the manner specified
%       by the author or licensor (but not in any way that suggests that
%       they endorse you or your use of the work).
%
%     - Noncommercial. You may not use this work for commercial purposes.
%
%     - Share Alike. If you alter, transform, or build upon this work,
%       you may distribute the resulting work only under the same or
%       similar license to this one.
%
%   Any of the above conditions can be waived if you get permission
%   from the copyright holder.  Nothing in this license impairs or
%   restricts the author's moral rights.
%
%   This software 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.

%computes the homogeneous transform matrix A of the pose vector v
function A=v2t(v)
    c=cos(v(3));
    s=sin(v(3));
    A=[c, -s, v(1) ;
        s,  c, v(2) ;
        0   0  1  ];
end

%computes the pose vector v from an homogeneous transform A
function v=t2v(A)
    v(1:2, 1)=A(1:2,3);
    v(3,1)=atan2(A(2,1),A(1,1));
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function f = cart2pol(f, varargin)
%CART2POL    Transforms a DISKFUN F(x,y) in Cartesian coordinates to a
%   CHEBFUN2 F(th, r) in polar coordinates, with -pi <= th <= pi, and 0 <=
%   r <=1. Working with F in this setting does not ensure that F is smooth
%   at the origin of the disk (i.e. r=0), and can reduce the
%   accuracy of some computations. It may be better to use the 'polar' flag
%   in DISKFUN, or to work directly with the CDR decomposition of F.
%
%   F = cart2pol(f, 'cdr') returns F(th, r) with -pi <= th <= pi and
%   -1 <= r <= 1. This is equivalent to the CDR decomposition of F.
%
% See also CDR.

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

% Use the CDR to build a chebfun2 with trigfuns for rows and chebfuns for
% columns

[c, d, r] = cdr(f);

% Restrict columns down to r in [0 , 1]:
if nargin < 2
    c = restrict(c, [0 1]);
end

% Form a chebfun2
f = c*d*r.';

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function preparedata

data = load('F:\Dropbox\Workspace\Projects2\Packages\DNN_Demo_2015_08_08\train_clean_multi_tiny.mat');

for i=1:length(data.noisy_list)
    [~,uttID] = fileparts(data.noisy_list{i});
    utt_list{i} = uttID;
end

for i=1:length(data.mfcc_clean_train)
    mfcc_clean_train{i} = single(data.mfcc_clean_train{i}(:,1:13)');
    mfcc_noisy_train{i} = single(data.mfcc_noisy_train{i}(:,1:13)');
end
for i=1:length(data.mfcc_clean_test)
    mfcc_clean_test{i} = single(data.mfcc_clean_test{i}(:,1:13)');
    mfcc_noisy_test{i} = single(data.mfcc_noisy_test{i}(:,1:13)');
end

[label_train phones] = alignment2phonePosterior(data.alignment_train, [], 3);
[label_test] = alignment2phonePosterior(data.alignment_test, [], 3);

save('mfcc_label.mat', 'utt_list', 'mfcc_clean_train', 'mfcc_noisy_train', 'mfcc_clean_test', 'mfcc_noisy_test', 'label_train', 'label_test', 'phones');

end

author: singaxiong

repo: SignalGraph


text:
y1=user1sym();
y2=user3sym();
y3=user5sym();
y4=user7sym();
y5=user10sym();
y6=user12sym();
%//SYMBOL BASED CONFIGURATION\\
u1=[y1,y2,y3,y4,y5,y6];
v1=[1,3,5,7,10,12];

%// CHIP BASED CONFIGURATION\\
u2=[s1,s2,s3,s4,s5,s6];
v2=[1,3,5,7,10,12];
plot(u1,v1,'r-+',u2,v2,'g-*');
grid on;
legend('u1,v1','u2,v2');
title('simulation for No of users vs BER in  SYMBOL BASED & CHIP BASED ');
 xlabel('NUMBER USERS');
 ylabel('BER');

author: Sable

repo: mcbench-benchmarks


text:
function [E,U,T] = acrobotEnergy(z,p)
% [E,U,T] = acrobotEnergy(z,p)
%
% This function computes the mechanical energy of the acrobot.
%
% INPUTS:
%   z = [4,n] = state vector
%   p = parameter struct:
%       .m1 = elbow mass
%       .m2 = wrist mass
%       .g = gravitational acceleration
%       .l1 = length shoulder to elbow
%       .l2 = length elbow to wrist
%
% OUTPUTS:
%   E = [1,n] = total mechanical energy
%   U = [1,n] = potential energy
%   T = [1,n] = kinetic energy
%
% NOTES:
%
%   states:
%       1 = q1 = first link angle
%       2 = q2 = second link angle
%       3 = dq1 = first link angular rate
%       4 = dq2 = second link angular rate
%
%   angles: measured from negative j axis with positive convention
%

q1 = z(1,:);
q2 = z(2,:);
dq1 = z(3,:);
dq2 = z(4,:);

[U,T] = autoGen_acrobotEnergy(q1,q2,dq1,dq2,p.m1,p.m2,p.g,p.l1,p.l2);

E = U+T;

end

author: MatthewPeterKelly

repo: OptimTraj


text:
% test for a simple anisotropic metric
name = 'fixed-2d';
name = 'fixed-3d';
name = 'varying-2d';

rep = 'results/anisotropic-fm/';
if not(exist(rep))
    mkdir(rep);
end

% create the main direction of the field
switch name
    case 'fixed-3d' % Fixed 3D tensor field
        n = 20;
        s = [n n n]; % size
        % main direction of the tensor
        u = [1 .1 .1];
        U = repmat( reshape(u,1,1,1,3),[n n n 1] );
    case 'fixed-2d' % Fixed 2D tensor field
        n = 40;
        s = [n n 1];
        % main direction of the tensor
        u = [1 .1];
        U = repmat( reshape(u,1,1,2),[n n 1] );
    case 'varying-2d' % spacially varying 2D field
        n = 200;
        % create 2D vector field
        s = [n n 1];
        U = randn(n,n,2);
        sigma = 30;
        for it=1:10
            U = perform_vf_normalization( perform_blurring(U, sigma) );
        end

end
U = perform_vf_normalization( U );

% test for various degree of anisotropy
aniso_list = [0.01 0.05 .1 .2 .5 1];

nstart = 1;
nstart = 8;
if nstart==1
    start_points = [n/2 n/2];
else
    start_points = rand(2,nstart);
    start_points(1,:) = rescale(start_points(1,:),.1,.9);
    start_points(2,:) = rescale(start_points(2,:),.1,.9);
    start_points = round( start_points*(n-1)+1 );
end
if strcmp(name(end-1:end), '3d')
    start_points(end/2) = ceil(s(end)/2);
end

for ianiso = 1:length(aniso_list)

    aniso = aniso_list(ianiso);

    % use cross product to compute the 2 remaining orthogonal directions
    if strcmp(name(end-1:end), '2d')
        % 3D field
        V = cat(3, -U(:,:,2), U(:,:,1)); % orthogonal vector
        T = perform_tensor_recomp(U,V, ones(n),ones(n)*aniso );
    else
        % 3D field
        U = cat(5, U, randn(s(1),s(2),s(3),3,2));
        U(:,:,:,:,3) = cross( U(:,:,:,:,1),U(:,:,:,:,2), 4 );
        U(:,:,:,:,3) = perform_vf_normalization( U(:,:,:,:,3) );
        U(:,:,:,:,2) = cross( U(:,:,:,:,1),U(:,:,:,:,3), 4 );
        U(:,:,:,:,2) = perform_vf_normalization( U(:,:,:,:,2) );
        Lambda = ones(s(1),s(2),s(3),3);
        Lambda(:,:,:,2:3) = aniso;
        T = perform_tensor_decomp_3d(U,Lambda);
    end

    D = perform_fast_marching(T, start_points);

    D1 = D(:,:,ceil(s(end)/2));
    if strcmp(name(end-1:end),'3d')
        T1 = T(:,:,ceil(s(end)/2), 1:2,1:2);
    else
        T1 = T;
    end

    clf;
    imageplot(D1);
    colormap jet(256);

    clf;
    hold on;
    options.sub = round(n/15);
    options.color = 'k';
    plot_tensor_field(T1, D1, options);
    h = plot(start_points(2,:),start_points(1,:), 'r.');
    set(h, 'MarkerSize', 20);
    colormap jet(256);

    saveas(gcf, [rep name '-aniso-' num2str(ianiso) '.png'], 'png');

end

author: alecjacobson

repo: gptoolbox


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [10000]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%
for i = 1:length(lambda_vec)
		lambda = lambda_vec(i)
		[theta] = trainLinearReg(X, y, lambda);
		error_train(i) = linearRegCostFunction(X, y, theta, 0);
		error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
end

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

end

author: scruel

repo: Notes-ML-AndrewNg


text:
function model = rmSearchFit_twoGaussianDoGBetaFixed(model,data,params,wProcess,t)
% rmSearchFit_twoGaussiansDoG - wrapper for 'fine' two DoG Gaussian fit
%
% model = rmSearchFit_twoGaussiansDoG(model,prediction,data,params);
%
% Second gaussian negative only.
%
% 2008/01 SOD: split of from rmSearchFit.

% now get original sigmas:
gridSigmas_unique = unique(params.analysis.sigmaMajor);
% add upper and lower limit:
expandRange    = params.analysis.fmins.expandRange;
gridSigmas = [0.001.*ones(expandRange,1); ...
    gridSigmas_unique; ...
    params.analysis.maxRF.*ones(expandRange,1)];

% fminsearch options
searchOptions = params.analysis.fmins.options;

vethresh = params.analysis.fmins.vethresh;

% convert to double just in case
params.analysis.X = double(params.analysis.X);
params.analysis.Y = double(params.analysis.Y);
params.analysis.allstimimages = double(params.analysis.allstimimages);

% amount of negative fits
nNegFit  = 0;
trends   = t.trends;
t_id     = t.dcid+1;

% initialize
if ~isfield(model,'rssPos')
    model.rsspos = zeros(size(model.rss));
end

if ~isfield(model,'rssNeg')
    model.rssneg = zeros(size(model.rss));
end

%-----------------------------------
% Go for each voxel
%-----------------------------------
progress = 0;tic;
for ii = 1:numel(wProcess),

    % progress monitor (10 dots)
    if floor(ii./numel(wProcess)*10)>progress,
        % print out estimated time left
        if progress==0,
            esttime = toc.*10;
            if floor(esttime./3600)>0,
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d hours.\n',...
                    mfilename,numel(wProcess),ceil(esttime./3600));
            else
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d minutes.\n',...
                    mfilename,numel(wProcess),ceil(esttime./60));
            end;
            fprintf(1,'[%s]:Nonlinear optimization (x,y,sigma):',mfilename);
        end;
        fprintf(1,'.');drawnow;
        progress = progress + 1;
    end;

    % volume index
    vi = wProcess(ii);
    vData = double(data(:,ii));

    % raw rss value (non-squared) - faster than sum(data(:,vi).^2)
    rawrss     = norm(vData);

    % reset tolFun: Precision of evaluation function.
    % We define RMS improvement relative to the initial raw 'no-fit' data
    % RMS. So, 1 means stop if there is less than 1% improvement on the fit:
    searchOptions.TolFun = params.analysis.fmins.options.TolFun.*rawrss;

    % start point from grid fit
    startParams = [model.x0(vi); ...
        model.y0(vi); ...
        model.s(vi);...
        model.s2(vi)];

    % tight search region [lowerbound upperbound]
    if params.analysis.scaleWithSigmas,
        step = params.analysis.relativeGridStep.*startParams(3);
        minstep = params.analysis.maxXY./2./params.analysis.minimumGridSampling;
        step = min(step,minstep);
        maxstep = params.analysis.maxXY./2./params.analysis.maximumGridSampling;
        step = max(step,maxstep);
    else
        step = params.analysis.maxXY./2./params.analysis.maximumGridSampling;
    end;
    boundary.xy    = startParams(1:2)*[1 1] + [-1 1;-1 1].*step.*expandRange;

    % gridSigmas==startParams(3), somehow this fails sometimes so we'll
    % look for the closest one. The min and max make sure the closestvalue
    % stays within the data range.
    [tmp,closestvalue] = sort(abs(gridSigmas-startParams(3)));
    closestvalue    = max(closestvalue(1),expandRange+1);
    closestvalue    = min(closestvalue,numel(gridSigmas)-(expandRange+1));
    boundary.sigma     = gridSigmas(closestvalue+[-1 1].*expandRange)';

    % boundary sigma two depends on boudary of sigma 1
    boundary.sigma2 = [boundary.sigma(1).*params.analysis.minSigmaRatio ...
                       params.analysis.sigmaRatioInfVal];

    % combine all boundaries
    bndParams2      = double([boundary.xy;...
        boundary.sigma;...
        boundary.sigma2]);

    % actual fitting routine
    if searchOptions.MaxIter>0
        outParams = ...
            fmincon(@(x) rmModelSearchFit_twoGaussiansDoGBetaFixed(x,vData,...
            params.analysis.X,...
            params.analysis.Y,...
            params.analysis.allstimimages,trends,params.analysis.betaRatioAlpha),...
            startParams,[],[],[],[],bndParams2(:,1),bndParams2(:,2),...
            @(x) distanceCon(x,startParams,step,params.analysis.minSigmaRatio),searchOptions);
    else
        outParams = startParams;
    end

    % make predictions
    Xv = params.analysis.X - outParams(1);   % positive x0 moves center right
    Yv = params.analysis.Y - outParams(2);   % positive y0 moves center up
    %betaRatio = params.analysis.betaRatioAlpha.*(outParams(3).^2./outParams(4).^2);
    betaRatio = (outParams(3)./outParams(4)).^params.analysis.betaRatioAlpha;
    rf(:,1) = exp( (Yv.*Yv + Xv.*Xv) ./ (-2.*(outParams(3).^2)) );
    rf(:,2) = exp( (Yv.*Yv + Xv.*Xv) ./ (-2.*(outParams(4).^2)) );
    rfNew = rf(:,1) - betaRatio.*rf(:,2);
    X = [params.analysis.allstimimages*rfNew trends];
    b    = pinv(X)*vData;
    % force positive fit (is a constraint in the
    % rmModelSearchFit_twoGaussianDoGBetaFixed, so it should return here to get the right values)
    b(1) = abs(b(1));

    % test 1 gaussian alone (must be done prior to with second gaussian)
    if searchOptions.MaxIter==0
        % without second gaussian
        % make only the positive and the negative rf
        % for the positive rf put all the negative rf-values to zero, for
        % the negative rf put all the positive rf-values to zero.

        rfBeta = rfNew.*b(1);
        posInd = rfBeta > 0;
        negInd = rfBeta < 0;
        rfPos = rfNew;
        rfNeg = rfNew;
        rfPos(negInd) = 0;
        rfNeg(posInd) = 0;
        XPos = [params.analysis.allstimimages*rfPos trends];
        XNeg = [params.analysis.allstimimages*rfNeg trends];
        rssPos = norm(vData-XPos*b).^2;
        rssNeg = norm(vData-XNeg*b).^2;
%         rfPos = max(rfBeta,0);
%         rfNeg = min(rfBeta,0);
%         XPos = [params.analysis.allstimimages*rfPos trends];
%         XNeg = [params.analysis.allstimimages*rfNeg trends];
%         rssPos = norm(vData-XPos).^2;
%         rssNeg = norm(vData-XNeg).^2;

    end

    % with second gaussian
    rss  = norm(vData-X*b).^2;

    % store results only if the first beta is positive, somehow fmincon
    % outputs negative fits. If the fit is negative keep old (grid) fit.
    if b(1)>0,
        model.x0(vi)   = outParams(1);
        model.y0(vi)   = outParams(2);
        model.s(vi)    = outParams(3);
        model.s_major(vi)    = outParams(3);
        model.s_minor(vi)    = outParams(3);
        model.s_theta(vi)    = 0;
        model.s2(vi)   = outParams(4);
        model.rss(vi)  = rss;
        model.b([1 t_id],vi) = b;
        if searchOptions.MaxIter==0
            model.rsspos(vi) = rssPos;
            model.rssneg(vi) = rssNeg;
        end
    else
        % change the percent variance explained to be just under the
        % current vethresh. So it counts as a 'coarse'-fit but can still be
        % included in later 'fine'-fits
        model.rss(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
        nNegFit = nNegFit + 1;
        if searchOptions.MaxIter==0
            model.rsspos(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
            model.rssneg(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
        end
    end;
end;

% end time monitor
et  = toc;
if floor(et/3600)>0,
    fprintf(1,'Done [%d hours].\n',ceil(et/3600));
else
    fprintf(1,'Done [%d minutes].\n',ceil(et/60));
end;
fprintf(1,'[%s]:Removed negative fits: %d (%.1f%%).\n',...
    mfilename,nNegFit,nNegFit./numel(wProcess).*100);

return;

%-----------------------------------
% make sure that the pRF can only be moved "step" away from original
% poisiton "startParams"
% For the two Gaussian model we add the additional constraint that the
% second Gaussian is at least twice as large as the first.
function [C, Ceq]=distanceCon(x,startParams,step,minRatio)
Ceq = [];
dist = x([1 2])-startParams([1 2]);
C(1) = sqrt(dist(1).^2+dist(2).^2) - step;
C(2) = minRatio - 0.001 - x(4)./x(3);
return;
%-----------------------------------

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function spline_test05 ( )

%*****************************************************************************80
%
%% TEST05 tests BASIS_MATRIX_BEZIER, BASIS_MATRIX_TMP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  n = 4;
  ndata = 4;
  nsample = 4;
  tdata = [ 0.0E+00, 0.0E+00, 1.0E+00, 1.0E+00 ]';
  ydata = [ 7.0E+00,  8.3333333E+00,   10.0E+00, 12.0E+00 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  BASIS_MATRIX_BEZIER sets up the basis\n' );
  fprintf ( 1, '    matrix for the uniform Bezier spline.\n' );

  mbasis = basis_matrix_bezier ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    TDATA, YDATA\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ndata
    fprintf ( 1, '%12f  %12f\n', tdata(i), ydata(i) );
  end

  left = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '    T, Spline(T)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ndata

    if ( i == 0 )
      tlo = tdata(1) - 0.5E+00 * ( tdata(2) - tdata(1) );
      thi = tdata(1);
    elseif ( i < ndata )
      tlo = tdata(i);
      thi = tdata(i+1);
    elseif ( ndata <= i )
      tlo = tdata(ndata);
      thi = tdata(ndata) + 0.5E+00 * ( tdata(ndata) - tdata(ndata-1) );
    end

    if ( i < ndata )
      jhi = nsample - 1;
    else
      jhi = nsample;
    end

    for j = 0 : jhi

      tval = ( ( nsample - j ) * tlo   ...
             + (           j ) * thi ) ...
             /   nsample;

      yval = basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval );

      if ( 0 < i & j == 0 )
        mark = '*';
      else
        mark = ' ';
      end

      fprintf ( 1, '%c  %12f  %12f\n', mark, tval, yval );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% eml_osem_test.m
% compare aspire and matlab versions of E-ML-OSEM
% Copyright 2001-07-23, Jeff Fessler, The University of Michigan

% generate data
if ~isvar('xfbp'), printm 'data'
	if has_aspire
		f.dir	= test_dir;
		f.wtf	= [f.dir 't,g.wtf'];
		f.wtr	= strrep(f.wtf, 'wtf', 'wtr');
		f.yi	= [f.dir 'yi.fld'];
		f.ci	= [f.dir 'ci.fld'];
		f.ri	= [f.dir 'ri.fld'];
		f.mask	= [f.dir 'mask.fld'];
	end
	em_test_setup; f.is3b = 0;
%	em3_test_setup; f.is3b = 1;
prompt
end

% block object for block iterative
if ~isvar('Gb'), printm 'Gb'
	f.nblock = 5;
	Gb = Gblock(G, f.nblock);
end

% matlab iterations
if ~isvar('xmat'), printm 'matlab E-ML-OSEM'
	f.niter = 8;
	f.pixmax = 6;
	f.precon = 'fast';
	f.precon = 'classic';

	if f.is3b % fix: needs work!
		shaper = @(proj) reshape(proj, n.n1*n.n2, n.a);
	else
		shaper = @(x) x;
	end
	xinit = ig.ones;
	[xmat precon] = eml_osem(xinit(ig.mask), Gb, ...
			shaper(yi), shaper(ci), shaper(ri), ...
			'niter', f.niter, 'pixmax', f.pixmax, ...
			'precon', f.precon);
	precon = ig.embed(precon);
	xmat = ig.embed(xmat);
	im plc 1 2
	im(1, xmat, 'Matlab E-ML-OSEM iterations'), cbar horiz
prompt
end

if ~has_aspire, return, end

% aspire iterations
if ~isvar('xasp'), printm 'aspire E-ML-OSEM'

	f.init	= [f.dir 'init.fld'];
	f.out	= [f.dir 'out.fld'];
	fld_write(f.init, xinit, 'check', 0)
	if exist(f.out, 'file'), delete(f.out), end

	if ~isvar('f.sys_type')
		f.sys_type = ['2z@' f.wtr '@-'];
	end
	f.saver	= 'stack,1';
	f.alg	= sprintf('osemc,%s,%d,%d,1', f.precon, f.nblock, sg.na);
	f.penal	= '-';
	f.method = sprintf('@%d@%s@%s', f.niter, f.alg, f.penal);
	f.com = sprintf(['i -chat 5 empl3 %s %s  %s %s 1 %s 1 %s %s' ...
			' %s %s 0 1 %g 0 -'], ...
		f.out, f.init, f.yi, f.ci, f.ri, f.sys_type, f.mask, ...
		f.method, f.saver, f.pixmax);
	os_run(f.com)

	xasp = double(fld_read(f.out));
	im(2, xasp, 'Aspire E-ML-OSEM iterations'), cbar horiz
end

if 0 % test aspire preconditioner
	pasp = squeeze(fld_read('p.fld'));
	im plc 1 2
	im(1, precon, 'precon'), cbar
	im(2, pasp, 'aspire'), cbar
	max_percent_diff(precon, pasp)
	minmax(precon-pasp)
return
end

if 1
	t = vcorrcoef(xasp, xmat);
	printf('corr. %g,%g', t, t-1)

	im plc 2 2
	im(1, xmat, 'xhat matlab'), cbar
	im(2, xasp, 'xhat aspire'), cbar
	im(3, (xasp-xmat)/max(xmat(:)), 'aspire-matlab'), cbar

	t1 = eql_obj(xmat, G, yi(:), ci(:), ri(:), [], ig.mask);
	t2 = eql_obj(xasp, G, yi(:), ci(:), ri(:), [], ig.mask);

	if im
		subplot(224)
		plot(0:f.niter, t1-t1(1), '-o', 0:f.niter, t2-t1(1), '-x')
		xlabel iteration, ylabel '\Phi change', legend('mat', 'asp', 4)
		title(sprintf('E-ML-OSEM, Nsubset=%d', f.nblock))
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [We, words] = loadWeWords(modelFile, modelFormat)
% modelFormat:
% 0 -- Matlab file
% 1 -- text file with a header line <numWords> <embDim>.  Subsequent lines has <word> <values>
% 2 -- text file with each line has <word> <values>
% 3 -- assume that there are two files modelFile.We, modelFile.words

  verbose=0;
  %% load from Matlab mat file
  if modelFormat==0
    %% load We
    load(modelFile, 'We', 'words');
    if ~exist('We', 'var')
      load(modelFile , 'allW');

      if isfield(allW, 'We')
        We = allW.We;
      else
        error('evaluateWordSim: no We in %s', modelFile);
      end
    end

    %% load vocab
    if ~exist('words', 'var')
      error('No words in the model file %s\n', modelFile);
    end
  elseif modelFormat==3 % modelFile.We, modelFile.words
    % words
    fid = fopen([modelFile '.words'], 'r');
    tmp = textscan(fid, '%s');
    words = tmp{1};
    fclose(fid);
    We = dlmread([modelFile '.We']);
    We = We';
  %% load from a text file
  else
    fid = fopen(modelFile, 'r');
    line = fgetl(fid);
    assert(ischar(line));
    if modelFormat==1 % header line <numWord> <embDim>
      [~, embDimStr] = strtok(line, ' ');
      embDim = str2double(embDimStr);
      line = fgetl(fid);
    else
      tokens = strsplit(line);
      embDim = length(tokens)-1;
    end

    if verbose==1
      fprintf(2, '# Reading file %s, embDim=%d\n', modelFile, embDim);
    end

    words = {};
    We = [];
    numWords = 0;
    while ischar(line)
      line = strtrim(line);
      if strcmp(line, '')
        line = fgetl(fid);
        continuel
      end

      [word, valueStr] = strtok(line, ' ');
      % word
      words = [words word];
      % emb
      We = [We sscanf(valueStr, '%f')];

      numWords = numWords+1;
      if verbose==1 && (mod(numWords, 10000)==0)
        fprintf(2, ' (%d)', numWords);
      end
      line = fgetl(fid);
    end

    if verbose==1
      fprintf(2, ' Done! Num words = %d\n', numWords);
    end
  end
end

author: lmthang


text:
function varargout = gallerysphere(name)
%CHEB.GALLERYSPHERE   Spherefun example functions.
%   F = CHEB.GALLERYSPHERE(NAME) returns a spherefun corresponding to NAME.
%   See the listing below for available names.
%
%   For example,  plot(cheb.gallerysphere('football')) plots the classic
%   icosahedral pattern of the Addias Telestar football (or soccer ball for
%   Americans).  For details of how each function is constructed try
%   'type +cheb/gallerysphere' or 'edit cheb.gallerysphere'.
%
%   [F,FA] = CHEB.GALLERYSPHERE(NAME) also returns the anonymous function
%   FA used to define the function. Some gallery functions are generated by
%   operations beyond the usual Spherefun constructor, so FA in those cases
%   is equal to F.
%
%   CHEB.GALLERYSPHERE with no input argument returns a function chosen at
%   random from the gallery.
%
%   CHEB.GALLERYSPHERE with no output argument creates a plot of the selected
%   function.
%
%   Available names:
%
%   deathstar   A function resembling the surface of the Death Star.
%   gaussian    Gaussian function on the sphere centered at Gauss's birth place.
%   geomag      Radial component of the International Geomagnetic Reference
%               field from the IGRF-12 model for 2015.
%   football    Icosahedral pattern found on a traditional soccer ball.
%   jet         A zonal jet stream over the mid-latitudes of the northern hemisphere.
%   moire       Moire pattern from waves generated at two point sources.
%   neamtu      Function created by Mike Neamtu for testing various spline
%               interpolation methods on the sphere (see Alfeld, Neamtu,
%               Schumaker, J. Comput. Appl. Math. 1996).
%   peaks       Peaks like function on the sphere taken from the geopeaks
%               function in the MATLAB mapping toolbox.
%   randn       Random linear combination of all real spherical harmonics
%               of exact degree 40.  The coefficients are generated from a
%               i.i.d. Gaussian (normal) distribution with std=1.
%   reprodkern  Reproducing kernel for spherical harmonics of degree 10
%               centered at (x,y,z) = (-1/sqrt(3),-1/sqrt(3),1/sqrt(3)).
%   soccerball  Same as football, but for the American users.
%   stripes     Alternating striped pattern.
%   vortices    Two antipodal vortices taken from  Nair, Cote, and
%               Staniforth (1999).
%   :)          A function to make you happy
%
%   Gallery functions are subject to change in future releases of Chebfun.
%
% See also CHEB.GALLERY, CHEB.GALLERYTRIG, CHEB.GALLERY2, CHEB.GALLERY3, CHEB.GALLERYDISK.

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

% If the user did not supply an input, return a function chosen at random
% from the gallery.
if ( nargin == 0 )
    names = {'football','jet','soccerball','deathstar', 'vortices'...
        'gaussian','reprodkern','geomag','peaks','neamtu','randn',...
        'moire','stripes',':)'};
    name = names{randi(length(names))};
end

type=1;                  % Default plotting type
addEarthPlot = 0;        % Flag on whether or not to plot the earth (1=yes)
viewAngle = [-37.5 30];  % Default viewing angle
clrmap = parula(64);     % Default colormap
clraxis = [];            % Color axis for the colormap; empty means use default

% The main switch statement.
switch lower(name)
    % The classic football (or soccerball) pattern
    case {'football','soccerball'}
        f = spherefun.sphharm(6,0) + sqrt(14/11)*spherefun.sphharm(6,5);
        fa = f;
        cntrlvl = -[0.25 0.25];
        type = 3;

    % A function resembling the surface of the Death Star
    case 'deathstar'
        fa = @(x,y,z) -(exp(-30*((y+sqrt(3)/2).^2 + x.^2 + (z-1/2).^2)) + exp(-100*z.^2));
        f = spherefun(fa);
        type = 1;
        clrmap = flipud(jet);
        viewAngle = [-35 8];

    % Two diametrically oposed vortices taken from
    % R. D. Nair, J. Cote, A. Staniforth, Cascade interpolation for
    % semi-Lagrangian advection over the sphere, Quart. J. Roy. Meteor.
    % Soc. 125 (1999) 1445-1468
    case 'vortices'
        rho = @(th) 3*(sin(th)); w = @(th) (3*sqrt(2)/2*sech(rho(th)).^2.*tanh(rho(th)))./(rho(th)+eps);
        fa = @(lam,th) -tanh(rho(th)/5.*sin(lam-20*w(th)));
        f = spherefun(fa);
        cntrlvl = [0 0];
        type = 3;
        clrmap = jet;

    % A simple Gaussian
    case 'gaussian'
        % Coordinates of Braunschweig, where Gauss was born;
        coords = [10.516667 52.266667]/180*pi;
        [xc,yc,zc] = sph2cart(coords(1),coords(2),1);
        fa = @(x,y,z,xc,yc,zc) exp(-20*((x-xc).^2 + (y-yc).^2 + (z-zc).^2));
        f = spherefun(@(x,y,z) fa(x,y,z,xc,yc,zc));
        addEarthPlot = 1;
        viewAngle = [80 15];

    % The reproducing kernel for all spherical harmonics of degree 20
    case 'reprodkern'
        [lam0,th0] = cart2sph(-1/sqrt(3),-1/sqrt(3),1/sqrt(3));
        f = spherefun(@(lam,th) sphRPK(lam,th,lam0,pi/2-th0,20));
        fa = f;
        viewAngle = [-50 5];

    % The IGRF-12 geomagnetic field
    case 'geomag'
        % Spherical harmonic coefficients for the IGRF-12 2015 Geomagnetic
        % Model; see http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
        c = [-29442.0, -1501.0, 4797.1, -2445.1, 3012.9, -2845.6, 1676.7,...
            -641.9, 1350.7, -2352.3, -115.3, 1225.6, 244.9, 582.0, -538.4,...
            907.6, 813.7, 283.3, 120.4, -188.7, -334.9, 180.9, 70.4, -329.5,...
            -232.6, 360.1, 47.3, 192.4, 197.0, -140.9, -119.3, -157.5, 16.0,...
            4.1, 100.2, 70.0, 67.7, -20.8, 72.7, 33.2, -129.9, 58.9, -28.9,...
            -66.7, 13.2, 7.3, -70.9, 62.6, 81.6, -76.1, -54.1, -6.8, -19.5,...
            51.8, 5.7, 15.0, 24.4, 9.4, 3.4, -2.8, -27.4, 6.8, -2.2, 24.2, 8.8,...
            10.1, -16.9, -18.3, -3.2, 13.3, -20.6, -14.6, 13.4, 16.2, 11.7,...
            5.7, -15.9, -9.1, -2.0, 2.1, 5.4, 8.8, -21.6, 3.1, 10.8, -3.3, 11.8,...
            0.7, -6.8, -13.3, -6.9, -0.1, 7.8, 8.7, 1.0, -9.1, -4.0, -10.5, 8.4,...
            -1.9, -6.3, 3.2, 0.1, -0.4, 0.5, 4.6, -0.5, 4.4, 1.8, -7.9, -0.7,...
            -0.6, 2.1, -4.2, 2.4, -2.8, -1.8, -1.2, -3.6, -8.7, 3.1, -1.5, -0.1,...
            -2.3, 2.0, 2.0, -0.7, -0.8, -1.1, 0.6, 0.8, -0.7, -0.2, 0.2, -2.2,...
            1.7, -1.4, -0.2, -2.5, 0.4, -2.0, 3.5, -2.4, -1.9, -0.2, -1.1, 0.4,...
            0.4, 1.2, 1.9, -0.8, -2.2, 0.9, 0.3, 0.1, 0.7, 0.5, -0.1, -0.3, 0.3,...
            -0.4, 0.2, 0.2, -0.9, -0.9, -0.1, 0.0, 0.7, 0.0, -0.9, -0.9, 0.4,...
            0.4, 0.5, 1.6, -0.5, -0.5, 1.0, -1.2, -0.2, -0.1, 0.8, 0.4, -0.1,...
            -0.1, 0.3, 0.4, 0.1, 0.5, 0.5, -0.3, -0.4, -0.4, -0.3, -0.8];
        % Compute the vertical component of the magnetic field
        f = 0*spherefun.sphharm(0,0);
        k = 1;
        for l=1:13
            f = f + -((l+1))*sqrt(4*pi/(2*l+1))*c(k)*spherefun.sphharm(l,0);
            k = k + 1;
            for m=1:l
                f = f - (-1)^m*((l+1))*sqrt(4*pi/(2*l+1))*...
                        (c(k)*spherefun.sphharm(l,m) + ...
                        c(k+1)*spherefun.sphharm(l,-m));
                k = k + 2;
            end
        end
        fa = f;
        type = 3;
        addEarthPlot = 1;
        cntrlvl = -60000:5000:60000;

    % A "peaks-like" function for the sphere
    case 'peaks'
        fa = @(x,y,z) 8*(1-x).^2.*exp(-4*(x - 0.059).^2 - 2*(y + 0.337).^2 - 2*(z + 0.940).^2) - ...
            30*(z/10 - x.^3 - y.^5) .* exp(-3*(x - 0.250).^2 - 2*(y - 0.433).^2 - 3*(z - 0.866).^2) + ...
            (20*y - 8*z.^3) .* exp(-2*(x + 0.696).^2 - 3*(y + 0.123).^2 - 2*(z - 0.707).^2) + ...
            (7*y - 10*x + 10*z.^3) .* exp(-3*(x - 0.296).^2 - 3*(y + 0.814).^2 - 3*(z + 0.5).^2);
        f = spherefun(fa);

    % A function developed by Mike Neamtu of Vanderbilt University to test
    % approximation schemes on the sphere.
    case 'neamtu'
        fa = @(x,y,z) 1 + x.^8 + exp(2*y.^3) + exp(2*z.^2) + 10*x.*y.*z;
        f = spherefun(fa);

    % A Guassian random function on the sphere suggested by Dmitry Belyaev
    % at Oxford.
    case 'randn'
        % Compute a random spherefun
        f = randnfunsphere(0.079,'monochromatic');
        fa = f;
        type = 4;
        clrmap = gray(2);
        % Set the color axis so that the transition from white to black is
        % exactly at the value zero.
        minmax = minandmax2est(f);
        clraxis = [-1 1]*max(abs(minmax));

    % An interesting Moire pattern generated by two sources on the sphere
    case 'moire'
        % Centers of the beacons
        boise = [-116.237651 43.613739]*pi/180;
        oxford = [-1.257778 51.751944]*pi/180;
        % ithaca = [-76.5 42.443333]*pi/180;
        % stellenbosh = [18.86 -33.92]*pi/180;
        [xb,yb,zb] = sph2cart(boise(1),boise(2),1);
        [xo,yo,zo] = sph2cart(oxford(1),oxford(2),1);
        % Pick the number of oscillations and make each of the "waves"
        % vanish at the anti-podal points from their centers.
        omega = besselroots(0,30); omega = omega(end)/2;
        % Use a combination of the J0 bessel functions centered at Boise
        % and Oxford to generate the Moire pattern.
        fa = @(x,y,z,omega) 2 + besselj(0,omega*sqrt((x-xb).^2+(y-yb).^2+(z-zb).^2)) + ...
            2 + besselj(0,omega*sqrt((x-xo).^2+(y-yo).^2+(z-zo).^2));
        f = spherefun(@(x,y,z) fa(x,y,z,omega));
        type = 1;
        addEarthPlot = 1;
        viewAngle = [32 8];

    % A zipper-like stripe pattern for the sphere
    case 'stripes'
        fa = @(x,y,z) (1 + cos(10*pi*x)).*exp(-exp(-20*z)) + (1 - cos(10*pi*x)).*exp(-exp(20*z));
        f = spherefun(fa);
        type = 3;
        cntrlvl = [1 1];
        viewAngle = [-36 8];
        clrmap = hot;

    % A zonal jet stream in the mid-latitudes of the northern hemisphere.
    case 'jet'
        fa = spherefun(@(lam,th) sin(th-(pi/4+0.01*cos(12*lam))).*exp(-300*(1-cos(th-(pi/4+0.01*cos(12*lam))))));
        f = spherefun(fa);
        type = 1;
        clrmap = jet;
        viewAngle = [32 8];
        addEarthPlot = 1;

    % A smiley face emoticon
    case ':)'
        fa = @(x,y,z) exp(-20*(4*(x+cos(pi/6)*sin(pi/4)).^2 + 4*(y+sin(pi/6)*sin(pi/4)).^2 + (z-cos(pi/4)).^2)) + ...
                      exp(-20*(4*(x+cos(-pi/6)*sin(pi/4)).^2 + 4*(y+sin(-pi/6)*sin(pi/4)).^2 + (z-cos(pi/4)).^2)) + ...
                   9*exp(-100*((x + 1).^2 - x.*z + (z - 1/2).^2 - 1/10).^2);
        f = spherefun(fa);
        viewAngle = [-75 20];

    otherwise
        error('CHEB:GALLERYSPHERE:unknown:unknownFunction', ...
            'Unknown function.')
end

% Only return something if there is an output argument.
if ( nargout > 0 )
    varargout = {f, fa};
    return;
end

ptitle = [name ', rank = ' num2str(length(f))];
% Determine the type of plot to make
if type==1
    % Plot the function with a grid
    surf(f,'grid','k-')
    axis off, title(ptitle)
elseif type==2
    contour(f)
    axis off, title(ptitle)
elseif type==3
    plot(f), hold on
    contour(f,cntrlvl,'k-'), hold off
    axis off, title(ptitle)
else
    plot(f)
    axis off, title(ptitle)
end

view(viewAngle);
colormap(clrmap);

% Set the color axis if required.
if ( ~isempty(clraxis) )
    caxis(clraxis);
end

if ( addEarthPlot )
    hold on, spherefun.plotEarth('w-'), hold off
end

end

function f = sphRPK(lam,th,lam0,th0,deg)
%SPHRPK  Reproducing kernel for the spherical harmonics of a given degree.
%   F = SPHRPK(LAM,TH,LAM0,TH0,DEG) is the reproducing kernel centered at
%   (LAM0,TH0) for the space of spherical harmonics of degree DEG.

t = cos(lam).*sin(th).*cos(lam0).*sin(th0) + sin(lam).*sin(th).*sin(lam0).*sin(th0) + cos(th).*cos(th0);
pl = legpoly(0:deg,[-1,1]);
[m,n] = size(t);
t = t(:);
c = ones(length(t),1)*((2*(0:deg)+1)/4/pi);
f = reshape(sum(c.*pl(t),2),m,n);

end


text:
function [val,lb,ub] = gridconstraints(fh,x_struct, p_struct, k, N, x, p, userdata)
gridConHandler = ocl.Constraint(userdata);
x = ocl.Variable.create(x_struct,x);
p = ocl.Variable.create(p_struct,p);

fh(gridConHandler,k,N,x,p);

val = gridConHandler.values;
lb = gridConHandler.lowerBounds;
ub = gridConHandler.upperBounds;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function p = bspeval(d,c,k,u)
% BSPEVAL  Evaluate B-Spline at parametric points
% -------------------------------------------------------------------------
% ADAPTATION of BSPEVAL from C Routine
% -------------------------------------------------------------------------
%
% Calling Sequence:
%
%   p = bspeval(d,c,k,u)
%
%    INPUT:
%
%       d - Degree of the B-Spline.
%       c - Control Points, matrix of size (dim,nc).
%       k - Knot sequence, row vector of size nk.
%       u - Parametric evaluation points, row vector of size nu.
%
%    OUTPUT:
%
%       p - Evaluated points, matrix of size (dim,nu)
%

nu = numel(u);
[mc,nc] = size(c);
                                                %   int bspeval(int d, double *c, int mc, int nc, double *k, int nk, double *u,int nu, double *p){
                                                %   int ierr = 0;
                                                %   int i, s, tmp1, row, col;
                                                %   double tmp2;
                                                %
                                                %   // Construct the control points
                                                %   double **ctrl = vec2mat(c,mc,nc);
                                                %
                                                %   // Contruct the evaluated points
p = zeros(mc,nu);                               %   double **pnt = vec2mat(p,mc,nu);
                                                %
                                                %   // space for the basis functions
N = zeros(d+1,1);                               %   double *N = (double*) mxMalloc((d+1)*sizeof(double));
                                                %
                                                %   // for each parametric point i
for col=1:nu                                    %   for (col = 0; col < nu; col++) {
                                                %     // find the span of u[col]
    s = findspan(nc-1, d, u(col), k);           %     s = findspan(nc-1, d, u[col], k);
    N = basisfun(s,u(col),d,k);                 %     basisfun(s, u[col], d, k, N);
                                                %
    tmp1 = s - d + 1;                           %     tmp1 = s - d;
    for row=1:mc                                %     for (row = 0; row < mc; row++)  {
        tmp2 = 0;                               %       tmp2 = 0.0;
        for i=0:d                               %       for (i = 0; i <= d; i++)
           tmp2 = tmp2 + N(i+1)*c(row,tmp1+i);  % 	tmp2 += N[i] * ctrl[tmp1+i][row];
        end                                     %
        p(row,col) = tmp2;                      %       pnt[col][row] = tmp2;
    end                                         %     }
end                                             %   }
                                                %
                                                %   mxFree(N);
                                                %   freevec2mat(pnt);
                                                %   freevec2mat(ctrl);
                                                %
                                                %   return ierr;
                                                %   }


text:
 function [err, sn, T1] = nufft1_err_mm(om, N1, J1, K1, type, alpha, beta)
%function [err, sn, T1] = nufft1_err_mm(om, N1, J1, K1, type, alpha, beta)
% Compute worst-case error for each input frequency for min-max 1D NUFFT.
% in:
%	om	[M,1]	digital frequency omega in radians
%	N1		signal length
%	J1		# of neighbors used per frequency location
%	K1		FFT size (should be > N1)
%	type		'sinc' 'diric' 'qr'
%	alpha	[L,1]	Fourier series coefficients of scaling factors
%			trick: or, "sn" if length N1
%	beta		scale gamma=2pi/K by this in Fourier series
%			typically is K/N (me) or 0.5 (Liu)
% out:
%	err	[M,1]	worst-case error over unit-norm signals
%	sn	[N,1]	scaling factors corresponding to alpha,beta
%	T1	[J,J]	T matrix
%
% Copyright 2001-12-7, Jeff Fessler, The University of Michigan

% if no arguments, give an example
if nargin < 4
	help(mfilename)
	N = 100; K = 2*N; gam = 2*pi/K;
	J = 14;
	om = gam * linspace(0,1,101);
	[alpha, beta, ok] = nufft_best_alpha(J, 2, K/N);
	if ~ok, alpha = []; beta = 0.5; end
	errd = nufft1_err_mm(om, N, J, K, 'diric');
	errs = nufft1_err_mm(om, N, J, K, 'sinc');
	errq = nufft1_err_mm(om, N, J, K, 'qr');
	semilogy(om/gam, errs, 'g-x', om/gam, errd, 'y-+', om/gam, errq, 'c-o')
	xlabel '\omega / \gamma', ylabel 'E_{max}(\omega)'
	legend('Tr sinc', 'Tr diric', 'QR approach')
return
end

if ~isvar('type') | isempty(type),	type = 'sinc'; end
if ~isvar('alpha') | isempty(alpha)
	alpha = [1];	% default Fourier series coefficients of scaling factors
end
if ~isvar('beta') | isempty(beta)
	beta = 0.5;	% default is Liu version for now
end

use_qr = logical(0);
if streq(type, 'sinc')
	use_true_diric = logical(0);
elseif streq(type, 'diric')
	use_true_diric = logical(1);
elseif streq(type, 'qr')
	use_qr = logical(1);
else
	error 'unknown type'
end

%
% see if 'best' alpha is desired
%
if ischar(alpha)
	if streq(alpha, 'uniform')
		alpha = [1];
		beta = 0.5;
	else
		if streq(alpha, 'best')
			L = 0;
		elseif streq(alpha, 'best,L=1')
			L = 1;
		elseif streq(alpha, 'best,L=2')
			L = 2;
		else
			error 'unknown alpha argument'
		end
		[alpha, beta, ok] = nufft_best_alpha(J1, L, K1/N1);
		if ~ok
			tmp = 'optimal alpha unknown for J=%d, K/N=%g, L=%d';
			warning(sprintf(tmp, J1, K1/N1, L))
			sn = ones(N1,1);
			err = nan;
			return
		end
	end
end

%
% if requested, return corresponding scaling factors too
%
if length(alpha) == N1	% trick: special way to give explicit sn's
	sn = alpha;
	if ~use_qr, error 'give sn only allowed for QR version', end
elseif nargout > 1 | use_qr
	sn = nufft_scale(N1, K1, alpha, beta);
end

%
% QR approach to error
%
if use_qr
	n = [0:N1-1]' - (N1-1)/2;
	[nn, jj] = ndgrid(n, 1:J1);
	gam1 = 2*pi/K1;
	C = exp(i * gam1 * nn .* jj) / sqrt(N1);
	S = spdiag(sn);
	A = S' * C;
	[Q,R] = qr(S' * C, 0);	% [N,J] compact QR decomposition

	do = col(om - gam1*nufft_offset(om, J1, K1));
	Db = exp(i * n * do') / sqrt(N1);	% [N,M]
	err = Db - Q * (Q' * Db);
	err = sqrt(sum(abs(err).^2,1))';	% [M]
return
end

tol = 0;
T1 = nufft_T(N1, J1, K1, tol, alpha, beta, use_true_diric);	% [J,J]
r1 = nufft_r(om, N1, J1, K1, alpha, beta, use_true_diric);	% [J,M]

%
% worst-case error at each frequency
%
Tr1 = T1 * r1;			% [J,M]
err = sum(conj(r1) .* Tr1).';	% [M,1]
err = min(real(err), 1);
err = sqrt(1 - err);		% caution: this "1 -" may cause numerical error

author: thomaskuestner


text:
clf; hold on
K = 10;
color = .7*[1 1 1];
dx = .9;
binedges = (0:5*K)+.5;
for ki = 1:K
    Data = poissrnd(2*ki,1,200);
    xpos = ki;
    p = boxyviolin(Data, binedges, xpos,dx, color)
    Data = randn(1,200)+3*K;
    xpos = ki;
    p = boxyviolin(Data, binedges, xpos,dx, [1 .8 .8])
end
shg

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Translate the coordinates back into the
% original image.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function B = getROIBBox(B, roi)

ratio = (roi([3 4]) - roi([1 2]))';
B = bsxfun(@times, B, ratio([1 2 1 2])');
B = bsxfun(@plus, B, roi([1 2 1 2]));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function T_hat = run_tc(params)
  T = params.T;
  Omega = params.Idx;

  %Omega = find(Idx);
  %Ak = T(Omega);

  %% ====================== Load data ==============================
  normalize              =        max(T(:))                     ;
  Xn                     =        T/normalize                   ;
  [n1,n2,n3]             =        size(Xn)                      ;

  %p                      =        0.3                           ;
  %Omega                  =        zeros(size(Xn))               ;
  %chosen                 =        randperm(n1*n2*n3,...
  %                                       round(p*n1*n2*n3))     ;
  %Omega(chosen)          =        1                             ;

  alpha                  =        1.05                             ;
  maxItr                 =        100                          ; % maximum iteration
  rho                    =        0.01                          ;

  A                      =        diag(sparse(double(Omega(:)))); % sampling operator
  b                      =        A * Xn(:)                     ; % available data
  bb                     =        reshape(b,[n1,n2,n3]);

  %% ================ main process of completion =======================
  [X] = LtSVD_TC(A ,b,rho,alpha ,[n1,n2,n3],maxItr, Xn(:), false, false);
  X                      =        X * normalize                 ;
  T_hat                  =        reshape(X,[n1,n2,n3])         ;

  %X_dif                  =        T-X                           ;
  %RSE                    =        norm(X_dif(:))/norm(T(:))     ;
end

author: andrewssobral


text:
clc;
close all;
clear all;

[num,txt,raw] = xlsread('data\final105.xls');

descriptions = raw(2:size(raw,1),2);
style_ratings = num(1:size(num,1),1);
comfort_ratings = num(1:size(num,1),4);
overal_ratings = num(1:size(num,1),5);

productids = cell2mat(raw(2:size(raw,1),7));
uniqueprodids =unique(productids);
outputdata = [];
outputdesc = {};

for productid = 1:size(uniqueprodids,1)

    rawdataindex = (find(productids==uniqueprodids(productid)));
    temp_style_rating = mean(style_ratings(rawdataindex));
    temp_comfort_ratings = mean(comfort_ratings(rawdataindex));
    temp_overal_ratings = mean(overal_ratings(rawdataindex));
    temp_individual_desc = descriptions(rawdataindex);
    temp_descriptions = strcat(temp_individual_desc{:},' ');
    outputdata = [outputdata; temp_style_rating, temp_comfort_ratings, temp_overal_ratings];
    outputdesc{end+1} = temp_descriptions;
end

hist(outputdata,500)
tic
%featurize_parallel(inputcellarray, nminFeatures, removeStopWords, doStem, grams, cores)
[featureVector,selectedheaderskeys] = featurize_parallel(outputdesc', 20, 0, 1, 1, 4);
toc
csvwrite('data\productsdesc.csv', featureVector);
csvwrite('data\productsvocab.csv', selectedheaderskeys);
csvwrite('data\productsratings.csv', outputdata);
%

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function cc_to_st_test01 ( )

%*****************************************************************************80
%
%% CC_TO_ST_TEST01 tests CC_TO_ST using a 1-based matrix.
%
%  Discussion:
%
%    This test uses a trivial matrix whose full representation is:
%
%          2  3  0  0  0
%          3  0  4  0  6
%      A = 0 -1 -3  2  0
%          0  0  1  0  0
%          0  4  2  0  1
%
%    The 1-based CC representation is
%
%      #  ICC  CCC  ACC
%     --  ---  ---  ---
%      1    1    1    2
%      2    2         3
%
%      3    1    3    3
%      4    3        -1
%      5    5         4
%
%      6    2    6    4
%      7    3        -3
%      8    4         1
%      9    5         2
%
%     10    3   10    2
%
%     11    2   11    6
%     12    5         1
%
%     13    *   13
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 July 2014
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 5;
  ncc = 12;

  acc = [ ...
    2.0,  3.0, ...
    3.0, -1.0,  4.0, ...
    4.0, -3.0,  1.0, 2.0, ...
    2.0, ...
    6.0, 1.0 ];
  ccc = [ ...
    1, 3, 6, 10, 11, 13 ];
  icc = [ ...
    1, 2, ...
    1, 3, 5, ...
    2, 3, 4, 5, ...
    3, ...
    2, 5 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_TO_ST_TEST01\n' );
  fprintf ( 1, '  Convert a 1-based CC matrix to ST format.\n' );
%
%  Print the CC matrix.
%
  cc_print ( m, n, ncc, icc, ccc, acc, '  The CC matrix:' );
%
%  Convert it.
%
  [ nst, ist, jst, ast ] = cc_to_st ( m, n, ncc, icc, ccc, acc );
%
%  Print the ST matrix.
%
  st_print ( m, n, nst, ist, jst, ast, '  The ST matrix:' );

  return
end


text:
function x = cvx_getdual( x )
x = x.dual_;

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function LC_map(lat,lon,maxlat,minlat,maxlon,minlon)

    %LC_MAP plot a map using a Lambert Conformal projection
    %
    %    LC_MAP(lat,lon,maxlat,minlat,maxlon,minlon)
    %
    %    Function to plot a map using a Lambert Conformal projection
    %    with two standard parallels which are chosen to be 1/4 of the
    %    vertical span of the map from the top (phi1) and bottom (phi2).
    %    The standard meridian is chosen to be the center meridian of
    %    the map.
    %
    %    where * lat & lon: array of latitudes and longitudes of map feature
    %	   * maxlat & minlat: maximum and minimum latitude of map
    %	   * maxlon & minlon: maximum and minimum longitude of map
    %	                      (remember: South & West are negative!)
    %
    %	The line type and line width can be set using the following
    %	global variables: "line_type" & "line_width".
    %	If these global variables are not set, it will use the
    %	following defaults: line_type = '-' & line_width = [0.5]

    report_this_filefun();

    global scale
    global phi0 lambda0 phi1 phi2
    global maxlatg minlatg maxlong minlong
    global line_type line_width

    ZG = ZmapGlobal.Data;
    % grab globals
    torad = ZG.torad;

    maxlatg = maxlat; minlatg = minlat;
    maxlong = maxlon; minlong = minlon;

    % set some constants
    scale = 1;

    % get the Standard Parallels and Center Coordinates
    phi1 = (minlat + ((maxlat - minlat) / 4)) * torad;
    phi2 = (maxlat - ((maxlat - minlat) / 4)) * torad;
    phi0 = (phi1 + phi2) / 2;
    lambda0 = ((minlon + maxlon) / 2) * torad;

    % find index of all points which are valid for this map & that are
    % not seperator points (ie: NaN)
    index = find(minlat < lat & lat < maxlat & minlon < lon & lon < maxlon &...
        isfinite(lat) & isfinite(lon));

    % convert all valid points to cartesian coordinates
    if index > 0
        [x(index) y(index)] = lc_tocart(lat(index),lon(index));
    end

    % reinsert the segment seperator points (ie: NaN)
    idxzero = find(x == 0 & y == 0);

    if idxzero > 0
        x(idxzero) = ones(size(x(idxzero))) * NaN;
        y(idxzero) = ones(size(y(idxzero))) * NaN;
    end

    % keep in memory if HOLD was on or off to put it back the way it was
    % after the plot
    if ishold
        hold_flag = 1;
    else
        hold_flag = 0;
    end

    lc_borde('-k',2)
    set(gca,'NextPlot','add')

    lc_grid(':',0.20)

    if isempty(line_type), line_type = '-k'; end
    if isempty(line_width), line_width = [0.5]; end
    plot(x,y,'-k','LineWidth',line_width)

    axis('equal')
    set(gca,'Visible','off')
    set(gcf,'PaperPosition',[1 .5 9 6.9545])

    % put HOLD back the way it was before this function was called
    if hold_flag
        set(gca,'NextPlot','add')
    else
        set(gca,'NextPlot','replace')
    end

end

function lc_borde(line_type,line_thick)

    %LC_PLOT_BORDER
    %
    %	LC_plot_border(line_type,line_thick)
    %
    %	Function to plot a border on the map generated by LC_MAP
    %	LINE_TYPE is any of the line types used by the PLOT command.
    %	LINE_THICK is the width of the line in "points".
    %	   (1 point = 1/72 inch).
    %	If no values are given, it will use the defaults: '-' & 2.
    %
    %	NOTE: The LC_MAP function has to have been called first
    %	to set some required global variables.

    global maxlatg minlatg maxlong minlong

    if nargin < 2
        line_thick = 2;
        if nargin < 1
            line_type = '-';
        end
    end

    bdlon(1:100) = minlong:(maxlong-minlong)/99:maxlong;
    bdlon(101:200) = ones(1,100) * maxlong;
    bdlon(201:300) = maxlong:-(maxlong-minlong)/99:minlong;
    bdlon(301:400) = ones(1,100) * minlong;

    bdlat(1:100) = ones(1,100) * maxlatg;
    bdlat(101:200) = maxlatg:-(maxlatg-minlatg)/99:minlatg;
    bdlat(201:300) = ones(1,100) * minlatg;
    bdlat(301:400) = minlatg:(maxlatg-minlatg)/99:maxlatg;

    [xbd, ybd] = lc_tocart(bdlat,bdlon);

    plot(xbd,ybd,line_type,'LineWidth',line_thick)

end

function lc_grid(line_type,line_width)

    %LC_PLOT_GRID
    %
    %	LC_plot_grid(line_type,line_width)
    %
    %	Function to plot a grid on the map generated by LC_MAP
    %	This function sets the size of the grid automatically according
    %	to the scale of the map.
    %	LINE_TYPE is the grid line style as available for the PLOT command
    %	(ie: [ - | -- | : | -. ]). The default is [ - ].
    %	LINE_WIDTH is the line thickness to be used for the grid. It has
    %	to be a value > [0.01].
    %
    %	NOTE: The LC_MAP function has to have been called before
    %	you can use this function as it needs some global variables
    %	to be set.

    global maxlatg minlatg maxlong minlong

    if nargin < 2
        line_width = [0.01];
        if nargin < 1
            line_type = '-';
        end
    end

    dlat = maxlatg - minlatg;
    dlon = maxlong - minlong;

    % figure out the grid increment for latitude and longitude
    if fix(dlat / 5) >= 1
        latinc = infix(dlat / 5);
    else
        if 2.5 < round(dlat)  && round(dlat) <= 5
            latinc = 1;
        elseif 1.5 < round(dlat)  &&  round(dlat) <= 2.5
            latinc = 1/2;
        else
            latinc = 1/4;
        end
    end

    if fix(dlon / 5) >= 1
        loninc = infix(dlon / 5);
    else
        if 2.5 < round(dlon)  && round(dlon) <= 5
            loninc = 1;
        elseif 1.5 < round(dlon)  &&  round(dlon) <= 2.5
            loninc = 1/2;
        else
            loninc = 1/4;
        end
    end

    % figure out at which latitude/longitude the grid starts
    % This is the best trick I found to do this in one step!
    ylatgr1 = ceil(minlatg/latinc)*latinc;
    xlongr1 = ceil(minlong/loninc)*loninc;

    % plot the parallel grid & labels
    k = 1;
    while ylatgr1 + (latinc * (k-1)) < maxlatg
        xlatgr = minlong:dlon/99:maxlong;
        ylatgr = ones(1,100) .* (ylatgr1 + (latinc * (k-1)));
        [x, y] = lc_tocart(ylatgr,xlatgr);
        plot(x,y,'LineWidth',line_width,'LineStyle',line_type)
        text(x(1),y(1),[num2str(ylatgr(1)) blanks(1)],'Units','data',...
            'HorizontalAlignment','right','VerticalAlignment','middle',...
            'FontWeight','bold','FontSize',12)
        k = k + 1;
    end

    % plot the meridian grid & labels
    k = 1;
    while xlongr1 + (loninc * (k-1)) < maxlong
        xlongr = ones(1,100) .* (xlongr1 + (loninc * (k-1)));
        ylongr = minlatg:dlat/99:maxlatg;
        [x, y] = lc_tocart(ylongr,xlongr);
        plot(x,y,'LineWidth',line_width,'LineStyle',line_type)
        text(x(1),y(1),[num2str(xlongr(1),5) blanks(1)],'Units','data',...
            'HorizontalAlignment','center','VerticalAlignment','top',...
            'FontWeight','bold','FontSize',12)

        k = k + 1;
    end

end

author: CelsoReyes


text:
function [LCpL,Q,sumLCpL,QE,Cy,M,Cp,Cq,Lq] = spm_eeg_assemble_priors(L,Qp,Qe,ploton,h)
% Predict sensor level impact of sources in Qp given sensor noise Qe
% FORMAT [LCpL,Q,sumLCpL,QE,Cy,M,Cp,Cq,Lq] = spm_eeg_assemble_priors(L,Qp,Qe,ploton,h)
% L       - lead fields
% Qp      - priors on source level dipole moment nAm/(mm2) per sample
% Qe      - sensor level variance in fT^2 per sample
% h       - optional hyperparameters that scale the variance components in
%           Qe and Qp (assume sensor followed by source level parameters)

% LCpL    - sensor level covariance components corresponding to the source
%           priors Qp
% Q       - sparse array over sources holding dipole moment density nAm/(mm2)
% sumLCpL - predicted sensor level variance due to sources (in fT^2)
% QE      - predicted sensor level noise variance (in fT^2)
% Cy      - total sensor noise covariance predicted : Cy = QE+sumLCpL
% M       - MAP estimator : M = Cp*L'*inv(Qe + L*Cp*L'))
% Cp      - total source covariance matrix
% Cq      - conditional source covariance - need to implement
% Lq      - cell array of the L*q (impact of each source component at sensor level)
%__________________________________________________________________________
% Copyright (C) 2010-2017 Wellcome Trust Centre for Neuroimaging

% Gareth Barnes
% $Id: spm_eeg_assemble_priors.m 7132 2017-07-10 16:22:58Z guillaume $

if nargin < 4
    ploton = ~spm('CmdLine');
end

if nargin < 5
    h=[];
end

Ne = length(Qe);
if iscell(Qp)
    Np = length(Qp); % number of source components
else
    Np = size(Qp,2);
end

if isempty(h)
    h = ones(Ne+Np,1);
end

if length(h)~=Ne+Np
    error('There must be as many hyperparameters as priors.');
end

he = h(1:Ne); % sensor level
hp = h(Ne + (1:Np)); % source level

Ns = size(L,2); % number sources

Q = sparse(Ns,Np);
Lq = {};

fprintf('Assembling %d prior components\n',Np);
for i = 1:Np
    if iscell(Qp)
        if isfield(Qp{i},'q')
            %Q(:,i) = Qp{i}.q; % patch amplitudes can be positive or negative
            q = Qp{i}.q;
            v = 1;
            if isfield(Qp{i},'v')
                v = sqrt(Qp{i}.v);
            end

            Q(:,i) = q*v;

            % Lq is the sensor level projection of the prior Q{i}.q
            Lq{i}.q = L*Q(:,i); %  supply an eigen mode in q
        else % no .q field, check it is 2D
            if size(Qp{i},1)~=size(Qp{i},2)
                disp('making Qp 2D');
                Qp{i} = diag(sparse(Qp{i}));
            end
        end
    end
end

%-Assemble empirical priors
%==========================================================================

LCpL  = {};

sumLCpL = zeros(size(L,1),size(L,1));

Cp  = sparse(0);
LCp = sparse(0);

for i = 1:Np

    if isfield(Qp{i},'q') % sparse diagonbal representation

        Q(:,i) = Q(:,i)*sqrt(hp(i)); % update moment by sqrt of variance scaling

        LQp = L*Q(:,i);
        LCpL{i} = LQp*LQp';
        Cp  = Cp + Q(:,i)*Q(:,i)'; %% Q is already scaled by root hp
        %LCp = LCp+LQp;
        LCp = LCp+L*Cp;
    else % full matrix version
        Qtmp = Qp{i}*hp(i);
        LCpL{i} = L*Qtmp*L';
        Cp  = Cp + Qtmp;
       % Q(:,i)=diag(Qtmp);
    end

    sumLCpL = sumLCpL + LCpL{i};

end

QE = zeros(size(Qe{1}));
for i=1:length(Qe)
    QE = QE + Qe{i}*he(i);
end

Cy = sumLCpL + QE;

% This is equivalent to M = Cp*UL'*inv(Qe + UL*Cp*UL'))
% with Cp the posterior source covariance (with optimal h values)
%Qsum=sum(Q(:,i).^2); % diagonal of posterior source cov matrix
%Cp=diag(Qsum);
M = Cp*L'/Cy;

% conditional variance (leading diagonal)
% Cq    = Cp - Cp*L'*iC*L*Cp;
%----------------------------------------------------------------------
%Cq    = Cp - sum(LCp.*M')'; in original

Cq = diag(Cp) - sum(LCp.*M')';

disp('end assemble');

if ploton
    subplot(3,1,1);
    imagesc(sumLCpL);colorbar;
    title('source prior');
    subplot(3,1,2);
    imagesc(QE);colorbar;
    title('sensor noise prior');
    subplot(3,1,3);
    imagesc(Cy);colorbar;
    title('Sensor covariance per sample (fT^2)');
end

author: spm


text:
function z = msign(x)
%% MSIGN modified sign function
%
% z = MSIGN(x) returns zero when abs(x)<sqrt(eps)

z = sign(x);
idx = (abs(x)< sqrt(eps));
z(idx) = 0;

repo: ifem


text:
function DEM_demo_double_well
% DEMO comparing DEM with particle filtering in the context of a bimodal
% conditional density.  This demonstrates a shortcoming of DEM in that it
% fails to represent the true density.

% get nonlinear state-space model
%==========================================================================
M      = spm_DEM_M('ssm');

% generate data (output)
%--------------------------------------------------------------------------
T      = 64;
U      = 8*sin(pi*(1:T)/16);
DEM    = spm_DEM_generate(M,U);

% EKF
%--------------------------------------------------------------------------
[kf_x] = spm_ekf(M,DEM.Y);

% PF
%--------------------------------------------------------------------------
[pf_x,P,Q,xQ] = spm_pf(M,DEM.Y);

% DEM
%--------------------------------------------------------------------------
DEM    = spm_DEM(DEM);
de_x   = DEM.qU.x{1};
tr_x   = DEM.pU.x{1};

spm_DEM_qU(DEM.qU,DEM.pU);

% Graphical comparison
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');

t      = 1:T;
subplot(2,2,1)
plot(t,pf_x,t,kf_x,':',t,de_x,t,tr_x)
legend({'PF','EKF','DEM','true'})
title('hidden state','FontSize',16)
axis([1 T -32 32])
axis square

subplot(2,2,2)
plot(t,tr_x - pf_x,t,tr_x - kf_x,':',t,tr_x - de_x)
legend({'PF','EKF','DEM'})
title('error','FontSize',16)
axis([1 T -32 32])
axis square

% Sample density
%--------------------------------------------------------------------------
subplot(2,2,3)
imagesc((1:T),xQ,Q)
axis xy square
title('sample density','FontSize',16)


text:
addpath(genpath('D:\GitHub\KiloSort2')) % path to kilosort folder
addpath('D:\GitHub\npy-matlab')

rootZ = 'H:\DATA\Spikes\Benchmark2015';

pathToYourConfigFile = 'D:\GitHub\KiloSort2\configFiles'; % take from Github folder and put it somewhere else (together with the master_file)
run(fullfile(pathToYourConfigFile, 'configFile384.m'))

ops.chanMap = 'H:\DATA\Spikes\Benchmark2015\forPRBimecToWhisper.mat';
rootH = 'H:\DATA\Spikes\temp\';

ops.NchanTOT = 129;
ops.trange      = [0 Inf]; % TIME RANGE IN SECONDS TO PROCESS

 % these settings overwrite any settings from config
ops.lam      = 10^2;   % weighting on the amplitude penalty (like in Kilosort1, but it has to be much larger)

ops.ThS      = [8 8];  % lower bound on acceptable single spike quality
ops.momentum = [20 400]; % number of samples to average over
ops.minFR    = 1/50; % minimum spike rate (Hz)

ops.sigmaMask  = 30;

ops.Nfilt       = 768; % max number of clusters
ops.nfullpasses = 1; % how many forward backward passes to do
ops.nPCs        = 3; % how many PCs to project the spikes into

ops.useRAM = 0;
ops.nNeighPC = 3;
ops.ccsplit = .975;

ops.NT                  = 64*1024+ ops.ntbuff;

dsets = {'20141202_all_es', '20150924_1_e', '20150601_all_s',...
    '20150924_1_GT', '20150601_all_GT', '20141202_all_GT'};

addpath(rootZ)

tic

for idk = [4 3 5] %:3 %1:length(dsets)
    fname = [dsets{idk} '.dat']; %'20141202_all_es.dat';

    savePath = fullfile(rootZ,  dsets{idk});

    ops.fbinary     = fullfile(savePath,  fname);
    ops.fproc       = fullfile(rootH, 'temp_wh2.dat'); % residual from RAM of preprocessed data

    % preprocess data
    rez = preprocessDataSub(ops);

    fname = fullfile(savePath, 'rez.mat');
    if exist(fname, 'file')
        dr = load(fname);
        rez.iorig = dr.rez.iorig;
        rez.ccb = dr.rez.ccb;
        rez.ccbsort = dr.rez.ccbsort;
    else
        rez = clusterSingleBatches(rez);
        save(fname, 'rez', '-v7.3');
    end

    rez = learnAndSolve8b(rez);

    rez = splitAllClusters(rez);

    fname = [dsets{idk} '.dat']; %'20141202_all_es.dat';
    [~, fn, ~] = fileparts(fname);

    isgood = rez.st3(:,4)<Inf;

    testClu = rez.st3(isgood,2);
    testRes = rez.st3(isgood,1);

    gtClu = LoadClu(fullfile(savePath, [fn '.clu.1']));
    fid   = fopen(fullfile(savePath, [fn '.res.1']), 'r');
    gtRes = int32(fscanf(fid, '%d'));
    fclose(fid);

    [allScores, allFPrates, allMissRates, allMerges] = ...
        compareClustering2(gtClu, gtRes, testClu, testRes, 1);

    rezToPhy(rez, savePath);

end

% save(fullfile(rootZ , 'allrez52718'), 'allrez', '-v7.3')

%%

%%

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function x = pwPoly2(tGrid,xGrid,t)
% x = pwPoly2(tGrid,xGrid,t)
%
% This function does piece-wise quadratic interpolation of a set of data.
%
% INPUTS:
%   tGrid = [1, 2*n-1] = time grid, knot idx = 1:2:end
%   xGrid = [m, 2*n-1] = function at each grid point in time
%   t = [1, k] = vector of query times (must be contained within tGrid)
%
% OUTPUTS:
%   x = [m, k] = function value at each query time
%
% NOTES:
%   If t is out of bounds, then all corresponding values for x are replaced
%   with NaN
%

nGrid = length(tGrid);
if mod(nGrid-1,2)~=0 || nGrid < 3
    error('The number of grid-points must be odd and at least 3');
end

% Figure out sizes
n = floor((length(tGrid)-1)/2);
m = size(xGrid,1);
k = length(t);
x = zeros(m, k);

% Figure out which segment each value of t should be on
edges = [-inf, tGrid(1:2:end), inf];
[~, bin] = histc(t,edges);

% Loop over each quadratic segment
for i=1:n
    idx = bin==(i+1);
    if sum(idx) > 0
        gridIdx = 2*(i-1) + [1,2,3];
        x(:,idx) = quadInterp(tGrid(gridIdx),xGrid(:,gridIdx),t(idx));
    end
end

% Replace any out-of-bounds queries with NaN
outOfBounds = bin==1 | bin==(n+2);
x(:,outOfBounds) = nan;

end

function x = quadInterp(tGrid,xGrid,t)
%
% This function computes the interpolant over a single interval
%
% INPUTS:
%   tGrid = [1, 3] = time grid
%   xGrid = [m, 3] = function grid
%   t = [1, p] = query times, spanned by tGrid
%
% OUTPUTS:
%   x = [m, p] = function at query times
%

% Rescale the query points to be on the domain [-1,1]
t = 2*(t-tGrid(1))/(tGrid(3)-tGrid(1)) - 1;

% Compute the coefficients:
a = 0.5*(xGrid(:,3) + xGrid(:,1)) - xGrid(:,2);
b = 0.5*(xGrid(:,3)-xGrid(:,1));
c = xGrid(:,2);

% Evaluate the polynomial for each dimension of the function:
p = length(t);
m = size(xGrid,1);
x = zeros(m,p);
tt = t.^2;
for i=1:m
    x(i,:) = a(i)*tt + b(i)*t + c(i);
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function r=findSubarrayWeights(T,d,a0)
%%FINDSUBARRAYWEIGHTS Given that the elements in an array (linear, planar)
%           have been grouped into subarrays and some type of element level
%           tapering might be applied, find the subarray-level tapering
%           that best approximates a tapering desired for alll elements in
%           an array. For example, one mgiht taper all the elements with a
%           Taylor tapering to get a good sum beam, but then want to design
%           subarray weights to be able to get a good difference beam.
%
%INPUTS: T A numSubarraysXnumEls matrix holding that breaks the elements
%          into subarrays and includes any element-level tapering.
%        d A numElsX1 vector of the desired element-level tapering. This
%          function provides a weighting for the subarrays to approximate
%          this.
%       a0 An optional numElX1 parameter. If one wishes to place a null in
%          the direction of a0 (often one might want to place a null in the
%          look direction when approximating a difference beam), then this
%          input can be provided. Otherwise, unconsitrained optimization is
%          performed. This parameter is the array manifold in the look
%          direction. For boresight, this will usually just be a numELX1
%          vector of ones. For an isotropic model, the vector is typically
%          a0=exp(-1j*2*pi*sum(bsxfun(@times,xyPoints,u),1)).' where
%          xyPoints is a 2XnumEl matrix of the locations of the elements in
%          the array.
%
%OUTPUTS: r The set of (possibly complex) weights for the subarrays to
%           approximate the given desired element-level tapering.
%
%This function implements the equations in [1]. The unconstrained
%optimization is just r=arg min_{r} norm(T*r-d)^2. The constrained
%optimization with a0 just adds the constraint that a0'*T*r'=0.
%
%Two approaches are given in the paper. The second approach considers
%emphasizing certain directions using a penalty function. However, when
%considering approximating the tapering for a Bayliss-weighted difference
%pattern, it was shown that this method was essentially the same as the
%simpler first method. Thus, this function just implements the first
%method.
%
%REFERENCES:
%[1] U. R. O. Nickel, "Subarray configurations for digital beamforming with
%    low sidelobes and adaptive interference suppression," in Record of the
%    IEEE International Radar Conference, Alexandria, VA, 8-11 May 1995,
%    pp. 714-719.
%
%August 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

T=T.';

numEls=size(T,1);

%Note that pinv(T)=inv(T'*T)*T'
TInv=pinv(T);

%If the constraint on placing a null in the look direction is imposed:
if(nargin>2&&~isempty(a0))
    r=TInv*(eye(numEls,numEls)-(a0*a0'*T*TInv)/(a0'*T*TInv*a0))*d;
else
    r=TInv*d;
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ output_data ] = deweightmem( input_data, weight_fun, oversample_rate, dim )
%DEWEIGHTMEM Make complex SAR uniformly weighted in one dimension
%
%    output_data = deweightmem(input_data, weight_fun, oversample_rate, dim)
%
%       Parameter name    Description
%
%       input_data        Array of complex values to deweight
%       weight_fun        Description of weighting applied.  Either a
%                            function handle to a function that takes a
%                            single argument (number of elements) and
%                            produces the weighting to apply, or a vector
%                            that is the weighting function sampled.
%       oversample_rate   Amount of sampling beyond the ImpRespBW in the
%                            processing dimension. (Default is Nyquist
%                            sampling = 1).
%       dim               Dimension over which to perform deweighting.
%                            Default is 1.
%       output_data       INPUT_DATA with normalization applied
%
% This implementation assumes that the data has already been "deskewed" and
% that the frequency support is centered.
%
% Author: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if ~exist('oversample_rate','var')
    oversample_rate = 1;
end
if ~exist('dim','var')
    dim = 1;
end

data_size = size(input_data,dim);
weight_size = round(data_size/oversample_rate); % Weighting only valid across ImpRespBW
if ~exist('weight_fun','var')||isempty(weight_fun) % No weighting passed in.  Do nothing.
    output_data = input_data;
    return;
elseif isa(weight_fun, 'function_handle')
    weighting = weight_fun(weight_size);
elseif isvector(weight_fun)
    weighting = interpft(weight_fun,weight_size);
end
% weighting = weighting * sqrt(mean(1./weighting.^2)); % We want total power maintained
weight_zp = ones(data_size,1); % Don't scale outside of ImpRespBW
weight_zp(floor((data_size-weight_size)/2)+(1:weight_size)) = weighting;
if dim==2
    weight_zp = weight_zp.';
end

% Divide out weighting in spatial frequency domain
output_data = fftshift(fft(input_data,[],dim),dim);
output_data = bsxfun(@rdivide,output_data,weight_zp);
output_data = ifft(ifftshift(output_data,dim),[],dim);

end

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

author: ngageoint


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

% Test script for illustration

% base scenario
t=0; T=1;
f = 0.03;
alpha = 0.075;
beta = 0.5;
nu = 0.1;
rho = 0;

sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);

lcoeff = 1; ucoeff = 1;
lowerbound = f-lcoeff*f;
upperbound = f + ucoeff*f;

xvals = lowerbound:.0005:upperbound;

zvals1d = density1d(xvals,sabrdensity,0,1);

% change rho
rho = -.9;
sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dm = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dm)
rho = .9;
sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dp = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dp)
figure('Color',[1 1 1]);
hold on;
plot(xvals,zvals1dm,'r','LineWidth',2); plot(xvals,zvals1d,'b','LineWidth',2); plot(xvals,zvals1dp,'g','LineWidth',2);
title('Plot of SABR densities for different \rho')
xlabel('0\leq x \leq .1')
ylabel('density')
legend('\rho=-0.9','\rho=0','\rho=0.9')
hold off;
rho = 0.0;

% change nu
nu = 0.05;
sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dm = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dm)
nu = .15;
sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dp = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dp)
figure('Color',[1 1 1]);
hold on;
plot(xvals,zvals1dm,'r','LineWidth',2); plot(xvals,zvals1d,'b','LineWidth',2); plot(xvals,zvals1dp,'g','LineWidth',2);
title('Plot of SABR densities for different \nu')
xlabel('0\leq x \leq .1')
ylabel('density')
legend('\nu=0.05','\nu=0.4','\nu=0.2')
hold off;
nu = 0.1;

% change beta
beta = 0.4;
sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dm = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dm)
beta = 0.8;
sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dp = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dp)
figure('Color',[1 1 1]);
hold on;
plot(xvals,zvals1dm,'r','LineWidth',2); plot(xvals,zvals1d,'b','LineWidth',2); plot(xvals,zvals1dp,'g','LineWidth',2);
title('Plot of SABR densities for different \beta')
xlabel('0\leq x \leq .1')
ylabel('density')
legend('\beta=0.4','\beta=0.5','\beta=0.6')
hold off;
beta = 0.5;

% change alpha
alpha = 0.05;
sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dm = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dm)
alpha = .1;
sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dp = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dp)
figure('Color',[1 1 1]);
hold on;
plot(xvals,zvals1dm,'r','LineWidth',2); plot(xvals,zvals1d,'b','LineWidth',2); plot(xvals,zvals1dp,'g','LineWidth',2);
title('Plot of SABR densities for different \alpha')
xlabel('0\leq x \leq .1')
ylabel('density')
legend('\alpha=0.05','\alpha=0.075','\alpha=0.1')
hold off;
alpha = 0.075;

author: Sable


text:
%FuNLMS
%Marko Stamenovic
%April 28, 2016

mus = [0.01 0.05 0.1 0.5 1 2];% 0.5 1];
for j = 1:length(mus)
    for i = 1:100
%% initialize values
% generate input signal
mu=0.13; %learning rate
M=128; %buffer size (num filter weights)
x=randn(10000,1); %input signal
x=x/max(x); %sample rate
fs=8000; %number of samples of the input signal
N=length(x); %length of input signal

% generate known filter coefficients
Pz=0.5*[0:127];  %linear coefficients
%Pz=randn(128,1); %random coefficients
ylim = max(Pz)*1.20;
ymin = min(Pz)-.2*max(Pz);
% generate filtered input signal == desired signal
d=conv(Pz,x); %input signal filtered by known filter Pz (primary path)

%% ESTIMATE SECONDARY PATH SIGNAL USING LMS %%
%generate dummy secondary path response Sz
Sz = Pz/2;
%run known signal through filter
yp=conv(Sz,x);
%initalize Sz hat values
Szh=zeros(M,1);

    for n=M:N
        ypvec=x(n:-1:n-M+1); %input has to be in reverse orxer has to be
        mu = 1/(ypvec'*ypvec);
        e(n)=yp(n)-Szh'*ypvec; %update error
        %plot(e)
        Szh=Szh+mu*ypvec*(e(n)); %update filter coefficient
    end
    Szh = abs(ifft(1./abs(fft(Szh))));

%% LMS FOR MAIN ANC %%
mu=2.6;
%initialize signal
y = zeros(N,1);
%filter input by learned filter to get x prime
xp = conv(Szh,x);
yp = conv(Szh,y);
%initialize delayed cancellation signal
ypvecLast = zeros(M,1);
yvecLast = zeros(M,1);
%initalize adaptive filter values
Az=zeros(M,1);
Bz=zeros(M,1);
%Make sure that x and d are column vectors
x=x(:);
d=d(:);
%LMS
    for n=M:N
        %flip buffer order
        yvec = y(n:-1:n-M+1);
        ypvec=yp(n:-1:n-M+1);
        xvec=x(n:-1:n-M+1)+.8*ypvec; %add feedback to input
        xpvec=xp(n:-1:n-M+1)+.8*ypvec; %add feedback to input
        %update mu
        mu = 1/(xvec'*xvec);
        %calculate output
        y(n) = Az'*xvec + Bz'*yvecLast;
        %update error
        e(n)=d(n)-y(n);
        %update adaptive filter weights
        Az=Az+mu*xpvec*(e(n));
        Bz=Bz+mu*ypvecLast*(e(n));
        %update delayed cancellation signal
        ypvecLast = ypvec;
        yvecLast = yvec;
        %draw the learned filter in realtime
%         plot(Pz)
%         hold on
%         plot(Az)
%         axis([0 inf ymin ylim])
%         title(sprintf('n=%f time=%fs error = %f mu = %f',n-M, (n-M)/fs, e(n),mu))
%         hold off
%         legend('Input coefficients','Learned Coefficients')
%         drawnow;
    end
    e=e(:);

%% PLOT RESULTS %%
figure
subplot(2,1,1)
plot(e)
title('Convergence Time in Cycles')
ylabel('Amplitude');
xlabel('Cycles');
legend('Error');
subplot(2,1,2)
stem(Pz)
hold on
stem(Az, 'r*')
title('Input Coefficients vs Learned Coefficients')
ylabel('Amplitude');
xlabel('Numbering of filter tap');
legend('Input Coefficients', 'learned coefficients')

author: markostam

repo: active-noise-cancellation


text:
function []=plot_uw(ifg_list);
%PLOT_UW plot some figures of intermediate unwrapping steps
%   PLOT_UW(IFG_LIST) - default is to plot the first
%
% Andy Hooper, March 2012
% modifications:
% DB    11/2014         Show the interferogram dates when plotted

if nargin<1
    ifg_list=1;
end

load uw_interp
load uw_grid
load ps2

drop_ifg = getparm('drop_ifg');
ifg_keep = [1:n_ifg]';
ifg_keep(drop_ifg)=[];
ifg_keep = ifg_keep(ifg_list);
if strcmpi(getparm('small_baseline_flag'),'y')
    ifg_str = [repmat('ifg ',length(ifg_keep),1)  num2str(ifg_keep) repmat(': ',length(ifg_keep),1) datestr(ifgday(ifg_keep,1)) repmat(' till ',length(ifg_keep),1)  datestr(ifgday(ifg_keep,2))];
else
    ifg_str = '';
end
ut=load('uw_space_time')

ph=ph(:,ifg_list);

n_ifg=size(ph,2);
ni=floor(sqrt(n_ifg+1));
nj=ceil(sqrt(n_ifg+1));
if ni*nj<n_ifg+1
    ni=ni+1;
end
[nrow,ncol]=size(Z);

figure
for i=1:n_ifg
     subplot(ni,nj,i)
     ifgw=ones(nrow,ncol)*(pi+pi/31.9);
     ifgw(nzix)=angle(ph(:,i));
     imagesc(ifgw,[-pi,pi+pi/31.9])
     axis off
     axis equal
     axis xy
     if ~isempty(ifg_str)
        title(ifg_str(i,:),'fontsize',15)
     end
end
colormap([jet(64);0,0,0]);
subplot(ni,nj,i+1)
imagesc(1,[-pi,pi])
colorbar('westoutside')
cla
axis off
set(gcf,'name','Resampled (Filtered) Phase')

figure
for i=1:n_ifg
    subplot(ni,nj,i)
    ifgw=angle(reshape(ph(Z,i),nrow,ncol));
    imagesc(ifgw,[-pi,pi])
    axis off
    axis equal
    axis xy
    if ~isempty(ifg_str)
        title(ifg_str(i,:),'fontsize',15)
     end
end
set(gcf,'name','Interpolated phase')
subplot(ni,nj,i+1)
imagesc(1,[-pi,pi])
colorbar('westoutside')
cla
axis off

load uw_space_time ifreq_ij jfreq_ij

if ~isempty(ifreq_ij)
    ifreq_ij=ifreq_ij(:,n_ifg);
    jfreq_ij=jfreq_ij(:,n_ifg);
    figure
    for i=1:n_ifg
        subplot(ni,nj,i)
        ifgw=zeros(nrow,ncol);
        clim=[min(ifreq_ij(:)),max(ifreq_ij(:))];
        clim(2)=clim(2)+diff(clim)/63.9;
        ifgw(nzix)=(ifreq_ij(:,i));
        ifgw(~nzix)=clim(2);
        imagesc(ifgw,clim);
        axis off
    %    axis equal
        axis xy
    end
    colormap([jet(64);0,0,0]);
    set(gcf,'name','North Phase Gradient')
    subplot(ni,nj,i+1)
    imagesc(1,clim)
    colorbar('westoutside')
    cla
    axis off
end

if ~isempty(jfreq_ij)
figure
for i=1:n_ifg
    subplot(ni,nj,i)
    ifgw=zeros(nrow,ncol);
    clim=[min(jfreq_ij(:)),max(jfreq_ij(:))];
    clim(2)=clim(2)+diff(clim)/63.9;
    ifgw(nzix)=(jfreq_ij(:,i));
    ifgw(~nzix)=clim(2);
    imagesc(ifgw,clim);
    axis off
%    axis equal
    axis xy
end
colormap([jet(64);0,0,0]);
set(gcf,'name','East Phase Gradient')
subplot(ni,nj,i+1)
imagesc(1,clim)
colorbar('westoutside')
cla
axis off
end

figure
ifgw=ones(nrow,ncol)*(pi+pi/31.9);
ifgw(nzix)=angle(ph(:,1));
imagesc([min(xy(:,2)),max(xy(:,2))],[min(xy(:,3)),max(xy(:,3))],ifgw,[-pi,pi+pi/31.9])
colormap([jet(64);0,0,0]);

plot_edges(edges,xy(:,2),xy(:,3),'w');
if isfield(ut,'shaky_ix')
    plot_edges(edges(ut.shaky_ix,:),xy(:,2),xy(:,3),'r');
end

axis off
axis xy

set(gcf,'name','Arcs')

end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function [nmps2] = mmps22nmps2(mmps2)
% Convert acceleration from millimeters per square-second to nanometers per
% second squared
% Chad A. Greene 2012
nmps2 = mmps2*1e+6;

repo: mcbench-benchmarks