text:
function [y, Dy, DDy, DDDy] = chebyshevInterpolate(f,t,d)
%
% This function uses barycentric interpolation to evaluate a chebyshev
% polynomial. Technical details are taken from the book: Approximation
% Theory and Approximation Practice by Lloyd Trefethen. The core equation
% that I use can be found as Theorem 5.2 in this book, and is called
% Barycentric Interpolation.
%
% Written by Matthew Kelly
% October 27, 2013
%   Updated: November 8, 2013
% Cornell University
%
% This function computes the value of the chebyshev polynomial that is
% defined by the vector f, at the inputs t. It can also be used to return
% the first and second derivatives of y with respect to t.
%
% INPUTS:
%   f = [KxN] value of the chebyshev polynomial at each of the chebyshev
%       points (these can be computed using chebyshevPoints.m). Each row
%       represents a single set of chebyshev node values for a single
%       state.
%   t = [1xM] vector of inputs to be evaluated (monotonically increasing)
%   d = [1x2] vector speficying the domain of the polynomial
%
% OUTPUTS:
%   y = [KxM] the value of the interpolant at each point in t
%   Dy = first derivative of y with respect to t
%   DDy = second derivative of y with respect to t
%   DDDy = third derivative of y with respect to t
%
% What is happening here, in plain english:
%
%   There are several Chebyshev points (or nodes) that are spread out
%   across the interval d using a special grid spacing. We will call these
%   points x. For each of the points in x, there is a corresponding value
%   of the chebyshev function, we'll call it f.
%
%   Let's say that we want to find the value of the approximation for some
%   input that is not in x. This can be done using a interpolation of the
%   values in f:
%                   y = c(t,x).*f
%
%   Note that the weighting terms are a function of both the desired input
%   and the chebyshev points. It turns out that there is a more stable and
%   efficient way to calculate this interpolation, which is roughly of the
%   form:
%                   y = (k(t,x).*f)/sum(k(t,x))
%
%   This code evaluates k(t,x) which it then uses to evaluate and return y.
%
%
% NOTE - The base algorithm (described above) is not defined for points in
% t that are close to the chebyshev nodes in x. If a point in t matches a
% point in x, then an additional claculation is required. If a point in t
% is very close to a gridpoint, then there is a slight loss of accuracy,
% which is more pronounced.
%

%Check to see if the t vector is on the proper domain:
if min(t) < d(1) || max(t) > d(2)
   disp('WARNING - the input to chebyshevInterpolate does not lie on the domain that was specified for the function')
   disp('    --> This can lead to unpredictable outputs and is not suggested.')
end

%Get the chebyshev points
[k,n] = size(f);
x = chebyshevPoints(n,d);
ONE1 = ones(k,1);
ONE2 = ones(1,length(t));

%Loop through each chebyshev node.
num = zeros(k,length(t));
den = zeros(k,length(t));
for i=1:n
    val = ONE1*(1./(t-x(i)));
    if mod(i,2)==1, val=-val; end;
    if i==1 || i==n
        num = num + 0.5*(f(:,i)*ONE2).*val;
        den = den + 0.5*(val);
    else
        num = num + (f(:,i)*ONE2).*val;
        den = den + val;
    end
end

%compute the solution:
y = num./den;

%Check for any values that were too close to nodes and correct them
nanIdx = isnan(y);
if sum(sum(nanIdx))>0
    nanRowIdx = max(nanIdx,[],1);
    y(:,nanRowIdx) = interp1(x',f',t(nanRowIdx)')';
end

%%%% Derivative Calculations %%%%
if nargout == 2
    Df = chebyshevDerivative(f,d);
    Dy = chebyshevInterpolate(Df,t,d);
elseif nargout == 3
    [Df, DDf] = chebyshevDerivative(f,d);
    Dy = chebyshevInterpolate(Df,t,d);
    DDy = chebyshevInterpolate(DDf,t,d);
elseif nargout == 4
    [Df, DDf, DDDf] = chebyshevDerivative(f,d);
    Dy = chebyshevInterpolate(Df,t,d);
    DDy = chebyshevInterpolate(DDf,t,d);
    DDDy = chebyshevInterpolate(DDDf,t,d);
end

end

% %  % This block of code provides a more efficient way of computing the
% %  % derivative of the function for very high order polynomial fits.
% %  %
% %  % It turns out that computing the differentiation matrix, which is done in
% %  % chebyshevDerivative, is an order n^2 calculation, which takes a long
% %  % time to compute for high order fits. This code is faster, but it is
% %  % not as accurate, especially near (or on) the gridpoints.
% %  %
% %  % To implement this code, put it between the lines:
% %  %    --> %%%% Derivative Calculations %%%%
% %  %    --> if nargout == 2
% %  %          CODE HERE
% %  %    --> end
% %  % Note that this would prevent the code from calculating higher order
% %  % derivatives.
% %
% %
% %
% %     %Loop through each chebyshev node again, slightly different math
% %     numDer = zeros(1,length(t));
% %     denDer = zeros(1,length(t));
% %     for i=1:n
% %         val = 1./(t-x(i)).^2;
% %         if mod(i,2)==1, val=-val; end;
% %         if i==1 || i==n
% %             numDer = numDer - 0.5*f(i).*(val);
% %             denDer = denDer - 0.5*(val);
% %         else
% %             numDer = numDer - f(i).*val;
% %             denDer = denDer - val;
% %         end
% %     end
% %
% %     %Compute the derivative of the approximation
% %     Dy = (numDer.*den - num.*denDer)./(den.^2);
% %
% %     %Check for any values that were to close to nodes
% %     %    -> NOTE - It is expensive to call chebyshevDerivative because it
% %     %       it computes a huge matrix just to take one or two derivatives.
% %     %    -> If speed is more important than accuracy, replace the contents
% %     %       of the folloging if statement with the code in comments after
% %     %       the end of the function. This uses a 4th order accurate central
% %     %       difference to compute the derivative.
% %     %
% %     nanIdx = isnan(Dy);
% %     nanSum = sum(nanIdx);
% %     if nanSum>0
% %
% %  %Use a 4th order central difference to get any values that were Nan
% %
% %         %This was determined experimentally:
% %         %  --> Too Small - run into machine precision
% %         %  --> Too Big - doesn't capture local shape well
% %         SMALL_NUMBER = 10^(-2*log10(length(f)));
% %
% %         %Get a vector of times that we need function values at
% %         timeDiff = [...
% %             t(nanIdx) - 2*SMALL_NUMBER;
% %             t(nanIdx) - SMALL_NUMBER;
% %             t(nanIdx) + SMALL_NUMBER;
% %             t(nanIdx) + 2*SMALL_NUMBER];
% %         timeDiffVec = reshape(timeDiff,4*nanSum,1)';
% %
% %         %Get the values of the function at each point in timeDiff
% %         funcDiffVec = chebyshevInterpolate(f, timeDiffVec,d)';
% %         funcDiff = reshape(funcDiffVec,4,nanSum);
% %
% %         %From Wikipedia: Finite difference coefficient
% %         %http://en.wikipedia.org/wiki/Finite_difference_coefficient
% %         coeff = [1/12, -2/3, 2/3, -1/12];
% %
% %         %Correct the Nan values
% %         Dy(nanIdx) = sum(coeff*funcDiff,1)/SMALL_NUMBER;
% %     end

author: MatthewPeterKelly

repo: dscTutorials


text:
%PREX_SOM   PRTools example on training SelfOrganizing Maps
%
% Show the training and plotting of 1- or 2-D Self-Organizing Maps.
%
help prex_som

delfigs
echo on
								% Set size of the SOM
	k = [5 1];
								% Set the number of iterations
   nrruns = [20 40 40];
								% Set desired learning rates
	eta = [0.5 0.1 0.1];
	                     % Set the neighborhood widths
   h = [0.6 0.2 0.01];
								% Generate one banana class:
	A = gendatb([100,100]);
	A = seldat(A,1);
								 % Train a 1D SOM:
  W = som(A,k);  % May take some time
								 % Show the results in a scatter plot
	figure(1); clf;
  scatterd(A); hold on;
	prplotsom(W);
	title('One-dimensional SOM');
	drawnow

								 % Train a 2D SOM:
	k = [5 5];
  W = som(A,k);  % Will take some time
								 % Show the results in a scatter plot
	figure(2); clf;
   scatterd(A); hold on;
	prplotsom(W);
	title('Two-dimensional SOM');
	drawnow

echo off
showfigs

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function [x,y,z] = genSquarePoints(app)
    method = app.DistributionDropDown.Value;
    switch method
        case 'Uniform'
            [x, y, z] = genUniformPoints(app);
        case 'Grid'
            [x, y, z] = genGriddedPoints(app);
    end
end

author: UMich-BipedLab


text:
function [retransformed_irf_estimates]=favar_retransX_irf_estimates(irf_estimates,transformationindex,levels)
% number of variables and shocks, and iterations (It-Bu)
[n1,n2]=size(irf_estimates);
ItBu=1;
% initialise
retransformed_irf_estimates=cell(size(irf_estimates));
% different treatment for different cases
if levels==1 % only cum sum, also for second differences
    for oo=1:n1 % variables
        for ll=1:n2 % shocks
            for ii=1:3 %for up.,median,low.
                if transformationindex(oo)==1 %1: no transformation, Level
                    retransformed_irf_estimates{oo,ll}(ii,:)=irf_estimates{oo,ll}(ii,:);
                elseif transformationindex(oo)==2
                    retransformed_irf_estimates{oo,ll}(ii,:)=cumsum(irf_estimates{oo,ll}(ii,:),2); %2: First Difference
%                     retransformed_irf_estimates{oo,ll}(ii,:)=irf_estimates{oo,ll}(ii,:);
                elseif transformationindex(oo)==3
                    %retransformed_irf_record{oo,ll}(ii,:)=cumsum(irf_record{oo,ll}(ii,:),2);    %3: Second Difference
                    retransformed_irf_estimates{oo,ll}(ii,:)=irf_estimates{oo,ll}(ii,:);
                elseif transformationindex(oo)==4
                    %retransformed_irf_record{oo,ll}(ii,:)=irf_record{oo,ll}(ii,:);
                    retransformed_irf_estimates{oo,ll}(ii,:)=exp(irf_estimates{oo,ll}(ii,:))-ones(ItBu,1);       %4: Log-Level
                elseif transformationindex(oo)==5
                    retransformed_irf_estimates{oo,ll}(ii,:)=cumsum(irf_estimates{oo,ll}(ii,:),2); %5: Log-First-Difference
%                     retransformed_irf_estimates{oo,ll}(ii,:)=exp(cumsum(irf_estimates{oo,ll}(ii,:),2))-ones(ItBu,1); %5: Log-First-Difference
                elseif transformationindex(oo)==6
                    %retransformed_irf_record{oo,ll}(ii,:)=cumsum(irf_record{oo,ll}(ii,:),2); %6: Log-Second-Difference
                    retransformed_irf_estimates{oo,ll}(ii,:)=exp(cumsum(irf_estimates{oo,ll}(ii,:),2))-ones(ItBu,1); %5: Log-First-Difference
                end
            end
        end
    end
elseif levels==2 % exp + cum sum
    for oo=1:n1 % variable
        for ll=1:n2 % shocks
            for ii=1:3 %for up.,median,low.
                if transformationindex(oo)==1 %1: no transformation, Level
                    retransformed_irf_estimates{oo,ll}(ii,:)=irf_estimates{oo,ll}(ii,:);
                elseif transformationindex(oo)==2
                    retransformed_irf_estimates{oo,ll}(ii,:)=cumsum(irf_estimates{oo,ll}(ii,:),2); %2: First Difference
                elseif transformationindex(oo)==3
                    retransformed_irf_estimates{oo,ll}(ii,:)=cumsum(irf_estimates{oo,ll}(ii,:),2);    % 6: Second Difference
                elseif transformationindex(oo)==4
                    retransformed_irf_estimates{oo,ll}(ii,:)=exp(irf_estimates{oo,ll}(ii,:))-ones(ItBu,1);       %4: Log-Level
                elseif transformationindex(oo)==5
                    retransformed_irf_estimates{oo,ll}(ii,:)=exp(cumsum(irf_estimates{oo,ll}(ii,:),2))-ones(ItBu,1); %5: Log-First-Difference
                elseif transformationindex(oo)==6
                    retransformed_irf_estimates{oo,ll}(ii,:)=exp(cumsum(irf_estimates{oo,ll}(ii,:),2))-ones(ItBu,1); %6: Log-Second-Difference
                end
            end
        end
    end
elseif levels==3 % exp + cum sum + cum sum for second differences types
    for oo=1:n1 % variables
        for ll=1:n2 % shocks
            for ii=1:3 %for up.,median,low.
                if transformationindex(oo)==1 %1: no transformation, Level
                    retransformed_irf_estimates{oo,ll}(ii,:)=irf_estimates{oo,ll}(ii,:);
                elseif transformationindex(oo)==2
                    retransformed_irf_estimates{oo,ll}(ii,:)=cumsum(irf_estimates{oo,ll}(ii,:),2); %2: First Difference
                elseif transformationindex(oo)==3
                    retransformed_irf_estimates{oo,ll}(ii,:)=cumsum(cumsum(irf_estimates{oo,ll}(ii,:),2),2);    % 6: Second Difference
                elseif transformationindex(oo)==4
                    retransformed_irf_estimates{oo,ll}(ii,:)=exp(irf_estimates{oo,ll}(ii,:))-ones(ItBu,1);       %4: Log-Level
                elseif transformationindex(oo)==5
                    retransformed_irf_estimates{oo,ll}(ii,:)=exp(cumsum(irf_estimates{oo,ll}(ii,:),2))-ones(ItBu,1); %5: Log-First-Difference
                elseif transformationindex(oo)==6
                    retransformed_irf_estimates{oo,ll}(ii,:)=exp(cumsum(cumsum(irf_estimates{oo,ll}(ii,:),2),2))-ones(ItBu,1); %6: Log-Second-Difference
                end
            end
        end
    end

else % do nothing, e.g. levels=0
    retransformed_irf_estimates=irf_estimates;
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [A,B] = hmxQRSVD(A,B,tol)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmxQRSVD.m                                    |
%|    #    |   VERSION    : 0.52                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.01.2019                                    |
%| ( === ) |   SYNOPSIS   : QR factorization and SVD recompression for    |
%|  `---'  |                low-rank matrices                             |
%+========================================================================+

if ~isempty(A)
    % QR Factorisation A = QA * RA;
    [QA,RA] = qr(A,0);

    % QR Factorisation Bt = QB * RB
    [QB,RB] = qr(B.',0);

    % Singular value decomposition U*S*V = RA * RB.'
    try
        [U,S,V] = svd(RA * RB.','econ');
    catch
        return
    end

    % Compression rank
    s = diag(S);
    n = sum( s./s(1) >= tol );

    % Recompression A = QA * U * S
    A = QA * (U(:,1:n) * S(1:n,1:n));

    % Recompression B = V' * QB^t
    B = V(:,1:n)' * QB.';
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [X,err,iter] = ksupport(A,B,k,opts)

% Solve the k support norm minimization problem by ADMM
%
% min_X 0.5*||vec(X)||_ksp^2, s.t. AX=B
% ---------------------------------------------
% Input:
%       A       -    d*na matrix
%       B       -    d*nb matrix
%       k       -    >0, integer, parameter
%       opts    -    Structure value in Matlab. The fields are
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       X       -    na*nb matrix
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 27/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

[d,na] = size(A);
[~,nb] = size(B);

X = zeros(na,nb);
Z = X;
Y1 = zeros(d,nb);
Y2 = X;

AtB = A'*B;
I = eye(na);
invAtAI = (A'*A+I)\I;

iter = 0;
for iter = 1 : max_iter
    Xk = X;
    Zk = Z;
    % update X
    temp = Z-Y2/mu;
    temp = prox_ksupport(temp(:),k,1/mu);
    X = reshape(temp,na,nb);
    % update Z
    Z = invAtAI*(-A'*Y1/mu+AtB+Y2/mu+X);
    dY1 = A*Z-B;
    dY2 = X-Z;
    chgX = max(max(abs(Xk-X)));
    chgZ = max(max(abs(Zk-Z)));
    chg = max([chgX chgZ max(abs(dY1(:))) max(abs(dY2(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', err=' num2str(err)]);
        end
    end

    if chg < tol
        break;
    end
    Y1 = Y1 + mu*dY1;
    Y2 = Y2 + mu*dY2;
    mu = min(rho*mu,max_mu);
end
err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);

author: canyilu

repo: LibADMM-toolbox


text:
function f  = replace_chromosome(intermediate_chromosome,pro,pop)
%% replace_chromosome(intermediate_chromosome,pro,pop)
% This function replaces the chromosomes based on rank and crowding
% distance. Initially until the population size is reached each front is
% added one by one until addition of a complete front which results in
% exceeding the population size. At this point the chromosomes in that
% front is added subsequently to the population based on crowding distance.

%
%  Copyright (c) 2009, Aravind Seshadri
%  All rights reserved.
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%  POSSIBILITY OF SUCH DAMAGE.

[N,V] = size(intermediate_chromosome);
switch pro
        case 1
        M = 2;
        V = 6;
    case 2
        M = 3;
        V = 12;
end

% Get the index for the population sort based on the rank
[temp,index] = sort(intermediate_chromosome(:,M + V + 1));

% Now sort the individuals based on the index
for i = 1 : N
    sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);
end

% Find the maximum rank in the current population
max_rank = max(intermediate_chromosome(:,M + V + 1));

% Start adding each front based on rank and crowing distance until the
% whole population is filled.
previous_index = 0;
for i = 1 : max_rank
    current_index = max(find(sorted_chromosome(:,M + V + 1) == i));
    if current_index > pop
        remaining = pop - previous_index;
        temp_pop = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        [temp_sort,temp_sort_index] = ...
            sort(temp_pop(:, M + V + 2),'descend');
        for j = 1 : remaining
            f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);
        end
        return;
    elseif current_index < pop
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
    else
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        return;
    end
    previous_index = current_index;
end

author: Sable

repo: mcbench-benchmarks


text:
function [VELperm, VELsystemC, VELsystemCT, rhsbcUx, rhsbcUy] = ...
      CurvedINSViscousSetUp2D(dt, nu, g0, BCfunction)

% function [VELperm, VELsystemC, VELsystemCT, rhsbcUx, rhsbcUy] = ...
%      CurvedINSViscousSetUp2D(dt, nu, g0, BCfunction)
% Purpose: build velocity system and boundary forcing terms

Globals2D;

% choose order to integrate exactly
Nint = ceil(3*N/2);

% build cubature nodes for all elements
CubatureOrder = 2*(Nint+1); cub = CubatureVolumeMesh2D(CubatureOrder);

% build Gauss node data for all element faces
NGauss = (Nint+1); gauss = GaussFaceMesh2D(NGauss);

% Convert boundary conditions to Dirichlet & Neumann
saveBCType = BCType;
ids = find(saveBCType==In | saveBCType==Wall | saveBCType==Cyl);  BCType(ids) = Dirichlet;
ids = find(saveBCType==Out);                    BCType(ids) = Neuman;

% Form inhomogeneous boundary term for rhs data (assumes time independent forcing)
[bcUx,bcUy,bcPR] = feval(BCfunction, gauss.x, gauss.y, gauss.nx, gauss.ny, gauss.mapI, gauss.mapO, gauss.mapW, gauss.mapC, 0, nu);

% Build pressure boundary condition forcing vector
[VELsystemBC] = CurvedPoissonIPDGbc2D();
rhsbcUx = VELsystemBC*bcUx(:); rhsbcUy = VELsystemBC*bcUy(:);

% Build velocity system
[VELsystem, mm] = CurvedPoissonIPDG2D();
VELsystem = g0*mm/(dt*nu) + VELsystem;

% Restore original boundary types
BCType = saveBCType;

% Find Reverse-Cuthill-Mckee ordering to reduce bandwith of velocity system
VELperm    = symrcm(VELsystem);

% Apply row and column permutation to velocity matrix
VELsystem  = VELsystem(VELperm,VELperm);

% Compute Cholesky factorization of velocity matrix
VELsystemC  = chol(VELsystem); VELsystemCT = transpose(VELsystemC);
return

repo: Aero-matlab


text:
% Test file for unbndfun/compose.m

function pass = test_compose(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

% Set the domain:
dom = [0 Inf];
domCheck = [0 1e2];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

%%%%%%%%%%%%%%%% Compose an UNBNDFUN with an operator (OP(F)) %%%%%%%%%%%%%%%%%%

opf = @(x) exp(-x);
opg = @(x) sin(exp(-x));
f = unbndfun(opf, struct('domain', dom));
g = compose(f, @sin);
gVals = feval(g, x);
gExact = opg(x);
err = gVals - gExact;
pass(1) = norm(err, inf) < 1e1*eps*get(g,'vscale');

%%%%%%%%%%%%%%%%%%% Compose two UNBNDFUNs (F + G) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opf = @(x) exp(-x);
opg = @(x) x.*exp(-x);
oph = @(x) (x+1).*exp(-x);
f = unbndfun(opf, struct('domain', dom));
g = unbndfun(opg, struct('domain', dom));
h = compose(f, @plus, g);
hVals = feval(h, x);
hExact = oph(x);
err = hVals - hExact;
pass(2) = norm(err, inf) < 1e1*eps*get(h,'vscale');

%%%%%%%%%%%%%%%%%%% Compose an UNBNDFUN with a BNDFUN (G(F)) %%%%%%%%%%%%%%%%%%%

opf = @(x) exp(-x);
opg = @(x) cos(x);
oph = @(x) cos(exp(-x));
f = unbndfun(opf, struct('domain', dom));
g = bndfun(opg, struct('domain', [-1 1]));
h = compose(f, g);
hVals = feval(h, x);
hExact = oph(x);
err = hVals - hExact;
pass(3) = norm(err, inf) < 1e1*eps*get(h,'vscale');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ w, t ] = circle_rule ( nt )

%*****************************************************************************80
%
%% CIRCLE_RULE computes a quadrature rule for the unit circle.
%
%  Discussion:
%
%    The unit circle is the region:
%
%      x * x + y * y = 1.
%
%    The integral I(f) is then approximated by
%
%      Q(f) = 2 * pi * sum ( 1 <= i <= NT ) W(i) * F ( cos(T(i)), sin(T(i)) ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NT, the number of angles to use.
%
%    Output, real W(NT), the weights for the rule.
%
%    Output, real T(NT), the angles for the rule.
%
  w(1:nt) = 1.0 / nt;
  for it = 1 : nt
    t(it) = 2.0 * pi * ( it - 1 ) / nt;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% fmri_glm_design_matrix
%
% Specify and/or create an fmri_glm_design_matrix object. This object is designed to contain
% all the information about a model necessary to create a design matrix
% and, when combined with an fmri_data object, to fit the model.
%
% The information contained in this object (its attributes) are designed to
% be very similar or identical to those defined by SPM in using
% spm_fmri_design.m. This should help with inter-operability and standards
% for storing information.
%
% See methods(fmri_glm_design_matrix) for things you can do with this.
% See the code of this function for additional details about all the
% fields.
%
% The general way to create/initialize CANlab objects is to call them with
% no or minimal input arguments. In this case, TR is required as the first
% input.  After that, pairs of arguments can be entered in the
%'fieldname', value format standard for Matlab inputs. If the fields are
% valid attributes of the object, they will be entered.
%
% For example:
%
% my_model = fmri_glm_design_matrix(2);
% creates an empty fmri_glm_design_matrix object with a TR of 2 (which is used to
% create a default b-spline basis set for the HRF.
%
% my_model = fmri_glm_design_matrix(2, 'nscan', [198 198 198]);
% creates an empty structure but assigns data to the field nscans, number
% of scans per session.
%
% my_model = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', 'onsets', ons_temperature, 'condition_names', names_temperature);
% does the above and also adds onsets and names to the appropriate place
% in the object structure.
% my_model.Sess(1).U(2).name : names for session 1, condition 2
% my_model.Sess(1).U(2).ons : onsets for session 1, condition 2 in UNITS(secs or TRs)
%
% add methods:
% There are some special methods to add specific types of data:
%
% 'onsets', {cell array of onsets for conditions within sessions}
%
% 'condition_names', {cell array of condition names} -- assumed to be the same for all sessions
%
% 'pm', '*name of condition to modulate*', '*modulator name*', {cell with modulators for each session}
%
% You can also specify different basis sets for different conditions, which
% SPM will not allow you to do.  Here is an example of a process of
% creating an fmri_glm_design_matrix object from onsets, etc., building it, and then
% replacing the basis set for one event type with another one.
%
% reportmod_model = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', 'onsets', ons_reportmod, ...
% 'condition_names', names_reportmod, 'pm', PM, 'pmnames', PM_names);
%
% reportmod_model = build(reportmod_model);
% plot(reportmod_model)
%
% % Generate a new basis set for the Anticipation conditions (condition 1)
% [xBF_hires, xBF] = fmri_spline_basis(2, 'length', 12, 'nbasis', 3, 'order', 3, 'plot');
% reportmod_model = replace_basis_set(reportmod_model, 1, xBF_hires);
% reportmod_model = build(reportmod_model);
% plot(reportmod_model)
%

%    SPM.Sess(s)
%            U: - Input structure array
%            C: - User specified covariate structure
%          row: - scan   indices for session s
%          col: - effect indices for session s
%           Fc: - F Contrast information for input-specific effects
%
%    SPM.xX
%            X: - design matrix
%           iH: - vector of H partition (indicator variables) indices
%           iC: - vector of C partition (covariates)          indices
%           iB: - vector of B partition (block effects)       indices
%           iG: - vector of G partition (nuisance variables)  indices
%         name: - cellstr of names for design matrix columns
%

%        3rd level
%        ------------------------------------------------------------------
%        SPM.Sess(s).U
%               dt: - time bin length {seconds}
%             name: - {1 x j} cell of names for each input or cause
%              ons: - (q x 1) onsets for q  trials {in UNITS}
%              dur: - (q x 1) durations for trials {in UNITS}
%                P: - Parameter stucture
%                u: - (t x j) inputs or stimulus function matrix
%              pst: - (1 x k) peristimulus times (seconds)
%
%
%        SPM.Sess(s).C
%
%                C: - [kx1 double] of user specified regressors
%             name: - {1xk} cellstr of regressor names
%
%
%        SPM.Sess(s).Fc
%
%                i: - F Contrast colums for input-specific effects
%             name: - F Contrast names  for input-specific effects
%
%
%            4th level
%            --------------------------------------------------------------
%            SPM.Sess(s).U(i).P(p)
%
%
%                 name: - parameter name
%                    P: - (q x 1) parameter matrix
%                    h: - order of polynomial expansion (0 = none)
%                    i: - sub-indices of U(i).u for plotting

classdef fmri_glm_design_matrix

    properties

        TR = NaN;

        % custom: not in SPM
        build_method
        history

        % 1st level
        % --------------------------------------------------------------------------
        % SPM.
        xY % : [1x1 struct] - data structure
        nscan % : [1xs double] - nscan(s) = number of scans in session s
        xBF % : [1x1 struct] - Basis function structure
        Sess % : [1xs struct] - Session structure array
        xX % : [1x1 struct] - Design matrix structure

    end % properties

    methods

        % Class constructor
        function obj = fmri_glm_design_matrix(varargin)
            %
            % [obj, cl_with_averages] = fmri_data(image_names, mask_image, varargin)
            %
            % Reads a set of image files and a mask image, and returns
            % an fmri_data object with data for all in-mask voxels.

            % ---------------------------------
            % Create empty fmri_data object, and return if no additional
            % arguments
            % ---------------------------------

%             if nargin == 0
%                 error('Must define TR, repetition time for scans, as first input.')
%             end

            if nargin > 0
                obj.TR = varargin{1};
            end

            obj.build_method = 'Separate sessions';

            obj.xY  = [];   % : [1x1 struct] - data structure
            obj.nscan = []; % : [1xs double] - nscan(s) = number of scans in session s
            obj.xBF = [];   % : [1x1 struct] - Basis function structure
            obj.Sess = [];  % : [1xs struct] - Session structure array
            obj.xX = [];

            % Default basis set
            %    ----------------------------------------------------------------------
            obj.xY.RT = obj.TR; % : - repetition time {seconds)

            obj.xBF = fmri_spline_basis(obj.TR, 0);
            % obj.xBF.name: - name of basis set
            % obj.xBF.length: - support of basis set {seconds}
            % obj.xBF.order: - order of basis set
            % obj.xBF.bf: - basis set matrix
            % obj.xBF.dt = TR ./ 16; % : - length of time bin {seconds}

            obj.xBF.T = 16; % : - number of time bins per scan
            obj.xBF.T0 = 1; % : - first time bin (see slice timing)
            obj.xBF.UNITS = [];  %: - 'scans'|'secs' - units in which onsets are specified
            obj.xBF.Volterra = 1; % : - 1|2 - order of [Volterra] convolution

            % Parameters for parametric modulation
            P = struct('name', '', 'P', [], 'h', 0, 'dur', [], 'i', []);

            % Onset structure
            U = struct('dt', obj.TR ./ obj.xBF.T , 'name', {}, 'ons', [], 'dur', 1, 'P', P, 'u', [], 'pst', []);

            % Covariate structure
            C = struct('C', [], 'name', {'User-specified regressors here'});

            % Session structure
            S = struct('U', U, 'C', C, 'row', [], 'col', [], 'Fc', []);

            obj.Sess = S;

            % Design structure
            obj.xX = struct('X', [], 'iH', [], 'iC', [], 'iB', [], 'iG', [], 'name', {});

            % The code below can be generic to any class definition
            % It parses 'fieldname', value pairs of inputs
            % and returns a warning if unexpected strings are found.

            if nargin == 1
                return
            end

            % all valid fieldnames
            valid_names = fieldnames(obj);

            for i = 1:length(varargin)
                if ischar(varargin{i})

                    % Look for a field (attribute) with the input name
                    wh = strmatch(varargin{i}, valid_names, 'exact');

                    % behaviors for valid fields
                    if ~isempty(wh)

                        obj.(varargin{i}) = varargin{i + 1};

                        % eliminate strings to prevent warnings on char
                        % inputs
                        if ischar(varargin{i + 1})
                            varargin{i + 1} = [];
                        end

                        % special methods for specific fields
                        switch varargin{i}

                        end

                    else
                        % Unique to this method: special subfields/other valid
                        % entries
                        % Try to add using the add method, which returns a
                        % warning if the field name is invalid.

                        obj = add(obj, varargin{i}, varargin{i + 1});

                        varargin{i + 1} = []; % eliminate to avoid confusing the parsing

                    end % not empty
                end % string input
            end % process inputs

        end % class constructor function

    end % properties

end % classdef

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [nuisance_ratio, rsquare_design, rsquare_nuis] = scn_component_rsquare(V, nuisanceX, designX)
% Print a table of r-square values (variance explained) for each of V data
% vectors by nuisance (mvmt, physio) and task-related predictors
%
% Designed to work with components
%
% :Examples:
% ::
%
%     % Typical operation
%     scn_component_rsquare(compscore, movement_params(1:157, :), X(1:157, :));
%
%     % No design
%     scn_component_rsquare(compscore, movement_params(1:157, :));
%
%     % Neither design nor nuisance, uses linear drift
%     scn_component_rsquare(compscore, []);
%
% ..
%    Tor Wager, Feb 2008
% ..

[t, m] = size(V);

if nargin < 3, designX = []; end

if isempty(nuisanceX)
    disp('Looking for nuisance covariates, but none found. Using linear drift.');
    nuisanceX = (1:t)';
end

disp('Variance in each component explained:')

% Nuisance
rsquare_nuis = get_rsquare(V, nuisanceX)';

if ~isempty(designX)
rsquare_design = get_rsquare(V, designX)';
else
    rsquare_design = .01 * ones(m, 1);
end

nuisance_ratio = rsquare_nuis ./ rsquare_design;

nuis_related = find(nuisance_ratio > 2);
task_related = find(nuisance_ratio < 1);

%rank_badness = sort(nuisance_ratio, 1, 'descend');
disp('All components');
print_matrix([(1:m)' rsquare_design rsquare_nuis nuisance_ratio], {'Comp.' 'R^2 Task' 'R^2 Nuisance' 'Ratio'});
fprintf('\n');

disp('Most task-related');
if isempty(designX)
    disp('NO DESIGN INFORMATION.');
else
print_matrix([task_related rsquare_design(task_related) rsquare_nuis(task_related) nuisance_ratio(task_related)], {'Comp.' 'R^2 Task' 'R^2 Nuisance' 'Ratio'});
end
fprintf('\n');

disp('Most nuisance-related');
print_matrix([nuis_related rsquare_design(nuis_related) rsquare_nuis(nuis_related) nuisance_ratio(nuis_related)], {'Comp.' 'R^2 Task' 'R^2 Nuisance' 'Ratio'});
fprintf('\n');

%fprintf('\tComp %3.0f : %3.0f%%\n', j, rsquare(j)*100);
end

function rsquare = get_rsquare(V, X)

    [t, m] = size(V);

    for j = 1:m
    % Center component to avoid counting intercept in r-square?
    % Don't have to, doesn't matter because var operator is 2nd moment
    y = V(:, j);
    b =  X \ y;

    fits = X * b;

    rsquare(j) = var(fits) / var(y);

    end

end


text:
function S=test(x)
% Edgar & Himmelblau, 1988
% x0 = [1, 2]'
% xo = [0, 0]'
% S(xo) = 0

S=4*x(1).^2-2.*x(1).*x(2)+x(2).^2;


text:
function imResult = blendMode_LinearDodge(A, B, offsetW, offsetH)
%% Linear Dodge blending mode: simply sums the values in the two layers.
%   Blending with white gives white. Blending with black does not change
%   the image.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_LinearDodge));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

C = blendMode_Add(A, B, offsetW, offsetH);

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end


text:
function [A,J2,A2] = create_adjacency_graph(X,method,param,force_sym)
%CREATE_ADJACENCY_GRAPH Summary of this function goes here
%   Detailed explanation goes here
if strcmp(method,'epsball')
    fh = @(Z) find_indices_by_epsilon(Z,param);
    [J2,A2] = create_adjacency_graph_internal(X,fh);
elseif strcmp(method,'nn')
    if 0 % my implementation
        fh = @(Z) find_indices_by_nn(Z,param);
        [J2,A2] = create_adjacency_graph_internal(X,fh);
    else % Matlab implementation is faster
        [D,I] = pdist2(X,X,'euclidean','Smallest',param+1);
        D(1,:) = [];
        I(1,:) = [];
        J2 = cell(size(D,2),1);
        A2 = cell(size(D,2),1);
        for i = 1:length(J2)
            J2{i} = I(:,i)';
            A2{i} = D(:,i)';
        end
    end
else
    disp('The method parameter can only take epsball or nn');
    A = []; J2 = []; A2 = [];
    return;
end
A = cellarr2sparse(J2,A2,length(J2),length(J2));
if force_sym
    A = max(A,A');
end

function [J2,A2] = create_adjacency_graph_internal(X,fh_find_indices)
step = 100;
n = size(X,1);
J2 = cell(n,1);
A2 = cell(n,1);
for i1 = 1:step:n
    i2 = i1 + step - 1;
    if (i2 > n)
        i2 = n;
    end;

    XX = X(i1:i2,:);
    D = calc_distance_matrix(XX,X);
    [Z,I] = sort(D,2);

    Z = Z(:,2:end); % remove itself as closest neighbor
    I = I(:,2:end);
    for i = i1:i2
        ind = fh_find_indices(Z(i-i1+1,:));
%         ind = find(Z(i-i1+1,:) < epsilon);
        jj = I(i-i1+1,ind);
        Z1 = Z(i-i1+1,ind);
        J2{i} = jj;
        A2{i} = Z1;
    end
end

function ind = find_indices_by_epsilon(Z,epsilon)
ind = find(Z < epsilon);

function ind = find_indices_by_nn(Z,nn)
ind = [1:min(nn,length(Z))];

author: zhouyuanzxcv

repo: Hyperspectral


text:
function spline_test16 ( )

%*****************************************************************************80
%
%% TEST16 tests SPLINE_CUBIC_SET, SPLINE_CUBIC_VAL2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  n = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST16\n' );
  fprintf ( 1, '  SPLINE_CUBIC_SET sets up a cubic spline;\n' );
  fprintf ( 1, '  SPLINE_CUBIC_VAL2 evaluates it.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Runge''s function, evenly spaced knots.\n' );

  for i = 1 : n

    t(i) = ( ( n - i     ) * (-1.0E+00)   ...
           + (     i - 1 ) * (+1.0E+00) ) ...
           / ( n     - 1 );

    y(i) =  frunge ( t(i) );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data to be interpolated:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of data values = %d\n', n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       T             Y\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '%14f  %14f\n', t(i), y(i) );
  end
%
%  Try boundary condition types 0, 1 and 2.
%
  for k = 0 : 2

    if ( k == 0 )

      ibcbeg = 0;
      ybcbeg = 0.0E+00;

      ibcend = 0;
      ybcend = 0.0E+00;

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Boundary condition 0 at both ends:\n' );
      fprintf ( 1, '  Spline is quadratic in boundary intervals.\n' );

    elseif ( k == 1 )

      ibcbeg = 1;
      ybcbeg = fprunge ( t(1) );

      ibcend = 1;
      ybcend = fprunge ( t(n) );

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Boundary condition 1 at both ends:\n' );
      fprintf ( 1, '  Y''(left) =  %f\n', ybcbeg );
      fprintf ( 1, '  Y''(right) = %f\n', ybcend );

    elseif ( k == 2 )

      ibcbeg = 2;
      ybcbeg = fpprunge ( t(1) );

      ibcend = 2;
      ybcend = fpprunge ( t(n) );

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Boundary condition 2 at both ends:\n' );
      fprintf ( 1, '  YP"(left) =  %f\n', ybcbeg );
      fprintf ( 1, '  YP"(right) = %f\n', ybcend );

    end

    ypp = spline_cubic_set ( n, t, y, ibcbeg, ybcbeg, ibcend, ybcend );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SPLINE"(T), F"(T):\n' );
    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '%14f  %14f\n', ypp(i), fpprunge(t(i)) );
    end

    left = 0;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  T, SPLINE(T), F(T), LEFT_IN, LEFT_OUT\n' );
    fprintf ( 1, '\n' );

    for i = 0 : n

      if ( i == 0 )
        jhi = 1;
      elseif ( i < n )
        jhi = 2;
      else
        jhi = 2;
      end

      for j = 1 : jhi

        if ( i == 0 )
          tval = t(1) - 1.0E+00;
        elseif ( i < n )
          tval = ( ( jhi - j + 1 ) * t(i)     ...
                 + (       j - 1 ) * t(i+1) ) ...
                 / ( jhi         );
        else
          if ( j == 1 )
            tval = t(n);
          else
            tval = t(n) + 1.0E+00;
          end
        end

        left_in = left;

        [ yval, ypval, yppval ] = ...
          spline_cubic_val2 ( n, t, y, ypp, left, tval );

        fprintf ( 1, '%14f  %14f  %14f  %6d  %6d\n', ...
          tval, yval, frunge ( tval ), left_in, left );

      end
    end

  end

  return
end


text:
function density = p02_density ( m, n, point )

%*****************************************************************************80
%
%% P02_DENSITY returns the density for problem 02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, real DENSITY(N), the mesh density at each point.
%
  density = ones ( 1, n );

  return
end


text:
function results = vl_test_aib(varargin)
% VL_TEST_AIB
vl_test_init ;

function s = setup()
s = [] ;

function test_basic(s)
Pcx = [.3 .3 0   0
       0   0   .2 .2] ;

% This results in the AIB tree
%
%  1 - \
%       5 - \
%  2 - /     \
%             - 7
%  3 - \     /
%       6 - /
%  4 - /
%
% coded by the map [5 5 6 6 7 1] (1 denotes the root).

[parents,cost] = vl_aib(Pcx) ;
vl_assert_equal(parents, [5 5 6 6 7 7 1]) ;
vl_assert_almost_equal(mi(Pcx)*[1 1 1], cost(1:3), 1e-3) ;

[cut,map,short] = vl_aibcut(parents,2) ;
vl_assert_equal(cut, [5 6]) ;
vl_assert_equal(map, [1 1 2 2 1 2 0]) ;
vl_assert_equal(short, [5 5 6 6 5 6 7]) ;

function test_cluster_null(s)
Pcx = [.5 .5   0   0
       0   0   0   0] ;

% This results in the AIB tree
%
%  1 - \
%       5
%  2 - /
%
%  3 x
%
%  4 x
%
% If ClusterNull is specified, the values 3 and 4
% which have zero probability are merged first
%
%  1 ----------\
%               7
%  2 ----- \   /
%           6-/
%  3 -\    /
%      5 -/
%  4 -/

parents1 = vl_aib(Pcx) ;
parents2 = vl_aib(Pcx,'ClusterNull') ;
vl_assert_equal(parents1, [5 5 0 0 1 0 0]) ;
vl_assert_equal(parents2(3), parents2(4)) ;

function x = mi(P)
% mutual information
P1 = sum(P,1) ;
P2 = sum(P,2) ;
x = sum(sum(P .* log(max(P,1e-10) ./ (P2*P1)))) ;

author: Cloud-CV

repo: object-proposals


text:
classdef testAmplificatorTensorNumericVsExplicitForSeqLam < testShowingError

    properties (Access = private)
        ampTensorNum
        ampTensorExp
        fiberDirection
        matValues
        C1
        C0
        Ch
        theta
    end

    properties (Access = protected)
       tol = 1e-10;
    end

    methods (Access = public)

        function obj = testAmplificatorTensorNumericVsExplicitForSeqLam()
            obj.createNumericalAmplificationTensor()
            obj.createExplicitAmplificationTensor()
        end

    end

    methods (Access = private)

        function createFiberDirection(obj)
           dir  = [1 0 0];
           obj.fiberDirection = Vector3D;
           obj.fiberDirection.setValue(dir);
           obj.fiberDirection.normalize();
        end

        function createNumericalAmplificationTensor(obj)
            obj.createFiberDirection()
            dir            = obj.fiberDirection;
            LevelOfFibers  = 4;
            FamilyName = 'HorizontalLaminate';
            LevelStr   = num2str(LevelOfFibers);
            name = strcat(FamilyName,LevelStr);
            printTopology  = true;
            iter           = 0;
            homogenizer    = NumericalFiberHomogenizer(dir,...
                             LevelOfFibers,name,...
                             printTopology,iter);
            obj.Ch         = homogenizer.getCh();
            obj.theta      = homogenizer.getVolume();
            P              = homogenizer.getAmplificatorTensor();
            obj.matValues  = homogenizer.getMaterialValues();
            obj.ampTensorNum = P;
        end

        function createExplicitAmplificationTensor(obj)
            obj.createMaterialTensors()
            c0 = obj.C0;
            c1 = obj.C1;
            dir = obj.fiberDirection;
            dirVal = dir.getValue;
            dirVal = [-dirVal(2),dirVal(1),0];
            dir.setValue(dirVal);
            t = obj.theta;
            Lam = AnisotropicLaminateHomogenizer(c0,c1,dir,t);
            P = Lam.getAmplificatorTensor();
            Ptens = CompliancePlaneStressTensor();
            Ptens.setValue(P);
            PtensV = Tensor2VoigtConverter.convert(Ptens);
            obj.ampTensorExp = PtensV;
        end

        function createMaterialTensors(obj)
            E1  = obj.matValues.E_plus;
            nu1 = obj.matValues.nu_plus;
            E0  = obj.matValues.E_minus;
            nu0 = obj.matValues.nu_minus;
            obj.C1  = IsotropicConstitutiveTensor(E1,nu1);
            obj.C0  = IsotropicConstitutiveTensor(E0,nu0);
        end

    end

    methods (Access = protected)

        function computeError(obj)
            va = obj.ampTensorNum.getValue();
            ta = obj.ampTensorExp.getValue();
            obj.error = norm(ta(:) - va(:))/norm(ta(:));
        end

    end

end


text:
%--- help for dsge/set_solution_to_companion ---
%
%  INTERNAL FUNCTION
%
%
%    Other functions named set_solution_to_companion
%
%       generic/set_solution_to_companion
%

author: jmaih


text:
function p05_story ( )

%*****************************************************************************80
%
%% P05_STORY prints the "story" for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This example is due to deBoor.\n' );
  fprintf ( 1, '  This data represents a temperature dependent property of titanium.\n' );

  return
end


text:
classdef DilationComputer < handle

    properties (Access = private)
        LHS
        RHS
    end

    properties (Access = private)
        mesh
        orientationVector
        dilation
    end

    methods (Access = public)

        function obj = DilationComputer(cParams)
            obj.init(cParams);
            obj.createDilationFun();
        end

        function rF = compute(obj)
            obj.computeLHS();
            obj.computeRHS();
            r = obj.solveSystem();
            s.mesh = obj.mesh;
            s.fValues = r;
            rF = P1Function(s);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh               = cParams.mesh;
            obj.orientationVector  = cParams.orientationVector;
        end

        function computeLHS(obj)
            K = obj.computeStiffnessMatrix();
            I = ones(size(K,1),1);
            obj.LHS = [K,I;I',0];
        end

        function K = computeStiffnessMatrix(obj)
            s.fun  = obj.dilation;
            s.mesh = obj.mesh;
            s.type = 'StiffnessMatrix';
            lhs = LHSintegrator.create(s);
            K = lhs.compute();
        end

        function createDilationFun(obj)
            obj.dilation = P1Function.create(obj.mesh, 1);
        end

        function computeRHS(obj)
            q = Quadrature.set(obj.mesh.type);
            q.computeQuadrature('CUBIC');
            gradT = obj.computeFieldTimesDivField(q);

            s.mesh = obj.mesh;
            s.type = 'ShapeDerivative';
            s.quadratureOrder = q.order;
            rhs  = RHSintegrator.create(s);
            rhsF = rhs.compute(gradT);
            rhsV = rhsF.fValues;
            obj.RHS = [rhsV;0];
        end

        function gradT = computeFieldTimesDivField(obj,q)
            a1    = obj.orientationVector{1};
            a2    = obj.orientationVector{2};
            aDa1  = a1.computeFieldTimesDivergence(q);
            aDa2  = a2.computeFieldTimesDivergence(q);
            s.quadrature = q;
            s.mesh       = obj.mesh;
            s.fValues    = -aDa1.fValues - aDa2.fValues;
            gradT = FGaussDiscontinuousFunction(s);
        end

        function u = solveSystem(obj)
            a.type = 'DIRECT';
            s = Solver.create(a);
            u = s.solve(obj.LHS,obj.RHS);
            u = u(1:end-1);
        end

    end

end


text:
function iclust = initialize_clusters(Ucell, Nk, type, Lx, Ly)

switch type
    case 'random'
        vs = randn(size(Ucell,1), Nk);
        vs = bsxfun(@rdivide, vs, sum(vs.^2,1).^.5 + 1e-8);% normalize activity vectors
        vs = single(vs);
        xs          = vs' * Ucell;
        [~, iclust] = max(xs,[],1);
    case 'Voronoi'

        xs = repmat(1:Lx, Ly, 1);
        ys = repmat((1:Ly)', 1, Lx);

        randx = rand(1, Nk) * Lx;
        randy = rand(1, Nk) * Ly;

        dx = repmat(xs(:), 1, Nk) - repmat(randx, numel(xs(:)), 1);
        dy = repmat(ys(:), 1, Nk) - repmat(randy, numel(ys(:)), 1);

        dxy = dx.^2 + dy.^2;
        [~, iclust] = min(dxy, [], 2);
    case 'squares'
        nsqrt = round(sqrt(Nk));

        xs = repmat(round(linspace(1, nsqrt, Lx)), Ly, 1);
        ys = repmat(round(linspace(1, nsqrt, Ly))', 1, Lx);
        iclust = xs + (ys-1) * nsqrt;

end

iclust = iclust(:)';

author: cortex-lab


text:
%MDL_QUADCOPTER Dynamic parameters for a quadrotor.
%
% MDL_QUADCOPTER is a script creates the workspace variable quad which
% describes the dynamic characterstics of a quadrotor flying robot.
%
% Properties::
%
% This is a structure with the following elements:
%
% nrotors   Number of rotors (1x1)
% J         Flyer rotational inertia matrix (3x3)
% h         Height of rotors above CoG (1x1)
% d         Length of flyer arms (1x1)
% nb        Number of blades per rotor (1x1)
% r         Rotor radius (1x1)
% c         Blade chord (1x1)
% e         Flapping hinge offset (1x1)
% Mb        Rotor blade mass (1x1)
% Mc        Estimated hub clamp mass (1x1)
% ec        Blade root clamp displacement (1x1)
% Ib        Rotor blade rotational inertia (1x1)
% Ic        Estimated root clamp inertia (1x1)
% mb        Static blade moment (1x1)
% Ir        Total rotor inertia (1x1)
% Ct        Non-dim. thrust coefficient (1x1)
% Cq        Non-dim. torque coefficient (1x1)
% sigma     Rotor solidity ratio (1x1)
% thetat    Blade tip angle (1x1)
% theta0    Blade root angle (1x1)
% theta1    Blade twist angle (1x1)
% theta75   3/4 blade angle (1x1)
% thetai    Blade ideal root approximation (1x1)
% a         Lift slope gradient (1x1)
% A         Rotor disc area (1x1)
% gamma     Lock number (1x1)
%
%
% Notes::
% - SI units are used.
%
% References::
% - Design, Construction and Control of a Large Quadrotor micro air vehicle.
%   P.Pounds, PhD thesis,
%   Australian National University, 2007.
%   http://www.eng.yale.edu/pep5/P_Pounds_Thesis_2008.pdf
% - This is a heavy lift quadrotor
%
% See also sl_quadrotor.

% 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

quadrotor.nrotors = 4;                %   4 rotors
quadrotor.g = 9.81;                   %   g       Gravity                             1x1
quadrotor.rho = 1.184;                %   rho     Density of air                      1x1
quadrotor.muv = 1.5e-5;               %   muv     Viscosity of air                    1x1

% Airframe
quadrotor.M = 4;                      %   M       Mass                                1x1
% Ixx = 0.082;
% Iyy = 0.082;
% Izz = 0.149;%0.160;
quadrotor.J = diag([0.082 0.082 0.149]);    %   I       Flyer rotational inertia matrix     3x3

quadrotor.h = -0.007;                 %   h       Height of rotors above CoG          1x1
quadrotor.d = 0.315;                  %   d       Length of flyer arms                1x1

%Rotor
quadrotor.nb = 2;                      %   b       Number of blades per rotor          1x1
quadrotor.r = 0.165;                  %   r       Rotor radius                        1x1

quadrotor.c = 0.018;                  %   c       Blade chord                         1x1

quadrotor.e = 0.0;                    %   e       Flapping hinge offset               1x1
quadrotor.Mb = 0.005;                 %   Mb      Rotor blade mass                    1x1
quadrotor.Mc = 0.010;                 %   Mc      Estimated hub clamp mass            1x1
quadrotor.ec = 0.004;                 %   ec      Blade root clamp displacement       1x1
quadrotor.Ib = quadrotor.Mb*(quadrotor.r-quadrotor.ec)^2/4 ;        %   Ib      Rotor blade rotational inertia      1x1
quadrotor.Ic = quadrotor.Mc*(quadrotor.ec)^2/4;           %   Ic      Estimated root clamp inertia        1x1
quadrotor.mb = quadrotor.g*(quadrotor.Mc*quadrotor.ec/2+quadrotor.Mb*quadrotor.r/2);    %   mb      Static blade moment                 1x1
quadrotor.Ir = quadrotor.nb*(quadrotor.Ib+quadrotor.Ic);             %   Ir      Total rotor inertia                 1x1

quadrotor.Ct = 0.0048;                %   Ct      Non-dim. thrust coefficient         1x1
quadrotor.Cq = quadrotor.Ct*sqrt(quadrotor.Ct/2);         %   Cq      Non-dim. torque coefficient         1x1

quadrotor.sigma = quadrotor.c*quadrotor.nb/(pi*quadrotor.r);         %   sigma   Rotor solidity ratio                1x1
quadrotor.thetat = 6.8*(pi/180);      %   thetat  Blade tip angle                     1x1
quadrotor.theta0 = 14.6*(pi/180);     %   theta0  Blade root angle                    1x1
quadrotor.theta1 = quadrotor.thetat - quadrotor.theta0;   %   theta1  Blade twist angle                   1x1
quadrotor.theta75 = quadrotor.theta0 + 0.75*quadrotor.theta1;%   theta76 3/4 blade angle                     1x1
quadrotor.thetai = quadrotor.thetat*(quadrotor.r/quadrotor.e);      %   thetai  Blade ideal root approximation      1x1
quadrotor.a = 5.5;                    %   a       Lift slope gradient                 1x1

% derived constants
quadrotor.A = pi*quadrotor.r^2;                 %   A       Rotor disc area                     1x1
quadrotor.gamma = quadrotor.rho*quadrotor.a*quadrotor.c*quadrotor.r^4/(quadrotor.Ib+quadrotor.Ic);%   gamma   Lock number                         1x1

quadrotor.b = quadrotor.Ct*quadrotor.rho*quadrotor.A*quadrotor.r^2; % T = b w^2
quadrotor.k = quadrotor.Cq*quadrotor.rho*quadrotor.A*quadrotor.r^3; % Q = k w^2

quadrotor.verbose = false;

author: petercorke

repo: robotics-toolbox-matlab


text:
function [h,hcb]=FakeColor(X,Y,Z,varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function h=FakeColor(X,Y,Z,clims,nc,colmap,plotcb,cbpos)
%
% Work-around for multiple colormaps
%
% ex.:
%     h=FakeColor(time,depth,u)
%     h=FakeColor(time,depth,u,0.6*[-1 1],10,'jet','right')
%     h=FakeColor(time,depth,u,0.6*[-1 1],10,'cool','right')
%
% INPUTS:
% X,Y,Z : Data to contour
% (Optional):
% clims: 1 X 2 vector with color limits [min max]
% nc: # color levels to contour (default 10)
% colmap: Colormap to use (ie 'cool'), default jet
% plotcb: true to plot, false to not plot colorbar
% cbpos: String specifiying position (see the colorbar documentation for
% options)
% Based on code at
% http://figuredesign.blogspot.com/2012/05/multiple-colormaps-in-matlab.html
% , posted by Sally Warner.
% Idea/method credit goes to Nick Beaird and Noel Pelland @ UW for figuring out
% how to do this. I just made it into a (hopefully) easy to use function.%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% AP - UW/APL - 10 Sept 2012
% apickering (at) apl (dot) washington (dot) edu
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Updates:
% 17 Spet 2012 - Make axes width same as if regular colorbar was inserted
% 20 June 2013 - Return handle to colorbar as well as axes
% 27 July 2013 - modified by Marcel Thielmann
%                   -added support for any
%                    colorbar location and changed the colorbar location
%                    identifier so that it complies with MATLAB standards
%                   -added the possibility to provide an own colormap as
%                    input
%                   -removed the axis ij command (can be done outside the
%                    function if desired)
%
% -COMMENT (Marcel Thielmann, 27 July 2013):
% When commands such as axis square are used afterwards, colorbar
% width/height does not fit axis widht/height anymore. I have not managed
% to automatically resize the colorbar once this is done, but the following
% workaround served the purpose (the plotboxpos function can be found on
% the file exchange:
% http://www.mathworks.com/matlabcentral/fileexchange/9615-plotboxpos )
%
% [c,cb] = FakeColor(x,y,z,CLIM,nc,'jet',1,'NorthOutside');
% axis square
% labels
% xlabel(c,'xlabel')
% ylabel(c,'ylabel')
%
% CPOS = plotboxpos(c);
% CBPOS = get(cb,'Position');

% set(cb,'Position',[CPOS(1) CBPOS(2) CPOS(3) CBPOS(4)])

%~~~~~~~~~~~~~~~~~
%%

nvars=numel(varargin);

% defaults if not specified
if nvars==0
    clims=[nanmin(Z(:)) nanmax(Z(:)) ];
    nc=15;
    colmap='jet';
    plotcb=1;
    cbpos='EastOutside';
elseif nvars==1
    clims=varargin{1};
    nc=15;
    colmap='jet';
    plotcb=1;
    cbpos='EastOutside';
elseif nvars==2
    clims=varargin{1};
    nc=varargin{2};
    colmap='jet';
    plotcb=1;
    cbpos='EastOutside';
elseif nvars==3
    clims=varargin{1};
    nc=varargin{2};
    colmap=varargin{3};
    plotcb=1;
    cbpos='EastOutside';
elseif nvars==4
    clims=varargin{1};
    nc=varargin{2};
    colmap=varargin{3};
    plotcb=varargin{4};
    cbpos='EastOutside';
elseif nvars==5
    clims=varargin{1};
    nc=varargin{2};
    colmap=varargin{3};
    plotcb=varargin{4};
    cbpos=varargin{5};
end

cvec=linspace(clims(1),clims(2),nc);

% make a "colormap"
%==============================================
if ischar(colmap)              % one of matlabs own builtin colormap given as string
    tcmap=eval([colmap '(nc)']);
else                           % own colormap given as matrix
    tcmap = zeros(nc,3);
    cvec_all = linspace(clims(1),clims(2),size(colmap,1));
    for i = 1:3
        tcmap(:,i) = interp1(cvec_all,colmap(:,i),cvec);
    end
end

% mimick a pcolor plot using contourf
[~,conhand]=contourf(X,Y,Z,cvec,'linecolor','none');hold on % ,'linecolor','none'
%
%get the patch object handles:
p=get(conhand,'children');
thechild=get(p,'CData');
cdat=cell2mat(thechild);
%
%loop through and manually set facecolor of each patch to the colormap you  made:
for i=1:length(cvec)
    set(p(cdat==cvec(i)),'Facecolor',tcmap(i,:))
end
%
h=gca; % get handle of axes

if plotcb
    %now you need to make a fake colorbar using the same trick:
    %define the colorbar axes location just above the subplot you are working with:

    % Make axes width same as if a regular colorbar was added; this way if
    % other axes have regular colorbars, all plot axes will (hopefully) line up

    cb = colorbar(cbpos);

    % get position of axes and colorbar for later use
    pos= get(gca,'position');
    CBPOS = get(cb,'position');
    cbprop = get(cb);
    % get axis locations and tickmarks for later use
    CBXTICK = get(cb,'XTick');
    CBYTICK = get(cb,'YTick');
    CBXAXLOC = get(cb,'XAxisLocation');
    CBYAXLOC = get(cb,'YAxisLocation');

    colorbar('off')
    set(gca,'position',pos);

    %contourf your fake colorbar:
    hcb             = axes('Position',CBPOS);

    % depending on the location, create a vertical or horizontal colorbar
    if ~isempty(strfind(cbpos,'North')) || ~isempty(strfind(cbpos,'South'))
        [~,conhand]     = contourf(cvec',[0 1],[cvec; cvec],cvec,'linecolor','none');
    else
         [~,conhand]     = contourf([0 1],cvec,[cvec' cvec'],cvec,'linecolor','none');
    end

    p               = get(conhand,'children');
    thechild        = get(p,'CData');
    cdat            = cell2mat(thechild);

    % AGAIN, loop through and manually set facecolor of each patch
    for i=1:length(cvec)
        set(p(cdat==cvec(i)),'Facecolor',tcmap(i,:))
    end
    set(hcb,'xtick',CBXTICK,'ytick',CBYTICK,'xaxislocation',CBXAXLOC,'yaxisLocation',CBYAXLOC)

else
    hcb = NaN;
end

%%


text:
function params = ReadParams(filename)

params = struct;
fid = fopen(filename,'r');
if fid < 0,
  error('Could not open parameter file %s for reading',filename);
end
while true,
  s = fgetl(fid);

  % end of file
  if ~ischar(s), break; end

  % comments
  if isempty(s) || ~isempty(regexp(s,'^\s*$','once')) || ...
        ~isempty(regexp(s,'^\s*#','once')),
    continue;
  end

  % split at first ,
  m = regexp(s,',','split','once');
  try
    name = m{1};
    if isempty(m{2}),
      val = {};
    else
      valcell = regexp(m{2},',','split');
      val = str2double(valcell);
      if any(isnan(val)),
        if numel(valcell) == 1,
          val = valcell{1};
        else
          val = valcell;
        end
      end
    end
    params.(name) = val;
  catch ME,
    warning('Unable to parse parameter line\n%s\nof %s:\n',s, ...
            filename,getReport(ME));
    continue;
  end

end

fclose(fid);

author: kristinbranson

repo: JAABA


text:
function node_num = fem2d_bvp_serene_node_num ( nx, ny )

%*****************************************************************************80
%
%% FEM2D_BVP_SERENE_NODE_NUM counts the number of nodes.
%
%  Discussion:
%
%    The program uses the finite element method, with piecewise serendipity
%    basis functions to solve a 2D boundary value problem over a rectangle.
%
%    The grid uses NX nodes in the X direction and NY nodes in the Y direction.
%
%    Both NX and NY must be odd.
%
%    Because of the peculiar shape of the serendipity elements, counting the
%    number of nodes and variables is a little tricky.  Here is a grid for
%    the case when NX = 7 and NY = 5, for which there are 29 nodes
%    and variables.
%
%     23 24 25 26 27 28 29
%     19    20    21    22
%     12 13 14 15 16 17 18
%      8     9    10    11
%      1  2  3  4  5  6  7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of X and Y grid values.
%    NX and NY must be odd and at least 3.
%
%    Output, integer NODE_NUM, the number of nodes and variables.
%
  node_num =   nx           * ( ny + 1 ) / 2 ...
           + ( nx + 1 ) / 2 * ( ny - 1 ) / 2;

  return
end


text:
%CREATECLIQUETREE Takes in a list of factors F, Evidence and returns a
%clique tree after calling ComputeInitialPotentials at the end.
%
%   P = CREATECLIQUETREE(F, Evidence) Takes a list of factors and creates a clique
%   tree. The value of the cliques should be initialized to
%   the initial potential.
%   It returns a clique tree that has the following fields:
%   - .edges: Contains indices of the nodes that have edges between them.
%   - .cliqueList: Contains the list of factors used to build the Clique
%   tree.
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function P = CreateCliqueTree(F, Evidence)

C.nodes = {};

V = unique([F(:).var]);

% Setting up the cardinality for the variables since we only get a list
% of factors.
C.card = zeros(1, length(V));
for i = 1 : length(V),

	 for j = 1 : length(F)
		  if (~isempty(find(F(j).var == i)))
				C.card(i) = F(j).card(find(F(j).var == i));
				break;
		  end
	 end
end

C.factorList = F;

% Setting up the adjacency matrix.
edges = zeros(length(V));

for i = 1:length(F)
	 for j = 1:length(F(i).var)
		  for k = 1:length(F(i).var)
				edges(F(i).var(j), F(i).var(k)) = 1;
		  end
	 end
end

cliquesConsidered = 0;

while cliquesConsidered < length(V)

	 % Using Min-Neighbors where you prefer to eliminate the variable that has
	 % the smallest number of edges connected to it.
	 % Everytime you enter the loop, you look at the state of the graph and
	 % pick the variable to be eliminated.

	 bestClique = 0;
	 bestScore = inf;
	 for i=1:size(edges,1)
		  score = sum(edges(i,:));
		  if score > 0 && score < bestScore
				bestScore = score;
				bestClique = i;
		  end
	 end

	 cliquesConsidered = cliquesConsidered + 1;
	 [F, C, edges] = EliminateVar(F, C, edges, bestClique);
end

% Pruning the tree.
C = PruneTree(C);

% We are incorporating the effect of evidence in our factor list.
for j = 1:length(Evidence),
	 if (Evidence(j) > 0),
		  C.factorList = ObserveEvidence(C.factorList, [j, Evidence(j)]);
	 end;
end;

% Assume that C now has correct cardinality, variables, nodes and edges.
% Here we make the function call to assign factors to cliques and compute the
% initial potentials for clusters.

P = ComputeInitialPotentials(C);

return

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function plot_lc(rho,eta,marker,ps,reg_param)
%PLOT_LC Plot the L-curve.
%
% plot_lc(rho,eta,marker,ps,reg_param)
%
% Plots the L-shaped curve of the solution norm
%    eta = || x ||      if   ps = 1
%    eta = || L x ||    if   ps = 2
% as a function of the residual norm rho = || A x - b ||.  If ps is
% not specified, the value ps = 1 is assumed.
%
% The text string marker is used as marker.  If marker is not
% specified, the marker '-' is used.
%
% If a fifth argument reg_param is present, holding the regularization
% parameters corresponding to rho and eta, then some points on the
% L-curve are identified by their corresponding parameter.

% Per Christian Hansen, IMM, 12/29/97.

% Set defaults.
if (nargin==2), marker = '-'; end  % Default marker.
if (nargin < 4), ps = 1; end       % Std. form is default.
np = 10;                           % Number of identified points.

% Initialization.
if (ps < 1 | ps > 2), error('Illegal value of ps'), end
n = length(rho); ni = round(n/np);

% Make plot.
loglog(rho(2:end-1),eta(2:end-1)), ax = axis;
if (max(eta)/min(eta) > 10 | max(rho)/min(rho) > 10)
  if (nargin < 5)
    loglog(rho,eta,marker), axis(ax)
  else
    loglog(rho,eta,marker,rho(ni:ni:n),eta(ni:ni:n),'x'), axis(ax)
    HoldState = ishold; hold on;
    for k = ni:ni:n
      text(rho(k),eta(k),num2str(reg_param(k)));
    end
    if (~HoldState), hold off; end
  end
else
  if (nargin < 5)
    plot(rho,eta,marker), axis(ax)
  else
    plot(rho,eta,marker,rho(ni:ni:n),eta(ni:ni:n),'x'), axis(ax)
    HoldState = ishold; hold on;
    for k = ni:ni:n
      text(rho(k),eta(k),num2str(reg_param(k)));
    end
    if (~HoldState), hold off; end
  end
end
xlabel('residual norm || A x - b ||_2')
if (ps==1)
  ylabel('solution norm || x ||_2')
else
  ylabel('solution semi-norm || L x ||_2')
end
title('L-curve')

repo: REST


text:
%% t_gifti_4
%

% Download from RDT
%{
 rdt = RdtClient('vistasoft');
 rdt.crp('/vistadata/anatomy/anatomyNIFTI');

 leftMeshFile = rdt.readArtifact('leftMesh',...
    'type','mat',...
    'destinationFolder',fullFolderName);
%}

%%  Set the vista data path and change into the gifti data directory

chdir(fullfile(vistaRootPath,'local'));
%% Create a gifti structure from a mrMesh mesh

% This loads a VISTASOFT msh structure
% First the original mesh
load(fullfile(vistaRootPath,'local','leftMesh.mat'));
g = gifti;
g.faces = int32(meshGet(msh,'triangles')' + 1);

% Choose either the initial unsmoothed vertices, or the smoothed vertices
% g.vertices = single(meshGet(msh,'initVertices')');
g.vertices = single(meshGet(msh,'vertices')');

% We don't know how this is used.
g.mat = eye(4,4);
% g.mat = rand(4,4);
% mrvNewGraphWin; h = plot(g);

% we take the green channel and use a gray scale map.  Their algorithm
% seems to normalize the color list and map it through the color map.
cdata = meshGet(msh,'colors');
cdata = cdata(2,:)';
gg.cdata = single(cdata);
mrvNewGraphWin; clf; colormap(gray); h = plot(g,gg);

%% Change the color map - I had some time on my hands.
colormap(cool); pause(0.5)
colormap(jet); pause(0.5)
colormap(redGreenCmap); pause(0.5)
colormap(blueyellowCmap); pause(0.5)

%% Change lighting and such
daspect([1,1,1]); view(45,30); axis tight
lightangle(45,30);
set(h,'SpecularColorReflectance',0,'SpecularExponent',50)
[az,el] = view;
g.mat = eye(4,4);
% g.mat = rand(4,4);
% mrvNewGraphWin; h = plot(g);

g.private.data{2}.metadata(1).name = 'AnatomicalStructurePrimary';
g.private.data{2}.metadata(1).value = 'CortexLeft';

% Also, DTI-Query will use the AnatomicalStructureSecondary metadata field
% to distinguish the surfaces that have been loaded.  While it's not
% necessary, it may be helpful to set this field as well:
% g.private.data{2}.metadata(2).name = 'AnatomicalStructureSecondary';
% g.private.data{2}.metadata(2).value = 'Pial'; % (or whatever you would like to tag the structure as)

% Then,
% save(g,'dhTest.gii');

% From Guillaume Flandin - storing metadata2
% g = gifti('file.gii');
% g.private.metadata(1).name  = 'AnatomicalStructurePrimary';
% g.private.metadata(1).value = 'CortexLeft';

%% Attach metadata to the gifti so Doug H. can read it. Name it as left.
% Helped by Guillaume Flandin email: gflandin@fil.ion.ucl.ac.uk

% First the left
load(fullfile(mrvDataRootPath,'anatomy','T1andMesh','Left_Mesh_Unsmoothed.mat'));
g = gifti;
g.faces = int32(meshGet(msh,'triangles')' + 1);

% Choose either the initial unsmoothed vertices, or the smoothed vertices
% g.vertices = single(meshGet(msh,'initVertices')');
g.vertices = single(meshGet(msh,'initVertices')');

% Set the flag to indicate this is a left hemisphere
g.private.data{2}.metadata(1).name = 'AnatomicalStructurePrimary';
g.private.data{2}.metadata(1).value = 'CortexLeft';
save(g,'left_gifti.gii');
% mrvNewGraphWin; h = plot(g); axis on; grid on

% Now the right
load(fullfile(mrvDataRootPath,'anatomy','T1andMesh','Right_Mesh_Unsmoothed.mat'));
g = gifti;
g.faces = int32(meshGet(msh,'triangles')' + 1);

% Choose either the initial unsmoothed vertices, or the smoothed vertices
% g.vertices = single(meshGet(msh,'initVertices')');
g.vertices = single(meshGet(msh,'initVertices')');

% Set the flag to indicate this is a left hemisphere
g.private.data{2}.metadata(1).name = 'AnatomicalStructurePrimary';
g.private.data{2}.metadata(1).value = 'CortexRight';

save(g,'right_gifti.gii');
% mrvNewGraphWin; h = plot(g);  axis on; grid on

% and actually, according to the specifications of the file format, you can
% store metadata for the complete file as above or for a given DataArray:
% g.private.data{2}.metadata(1).name = 'AnatomicalStructurePrimary';
% g.private.data{2}.metadata(1).value = 'CortexLeft';
%
% We need to figure out how to save these fields to work with dtiQuery.

%% End

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function display(anno)
fprintf(1,'\n %s = \n', inputname(1));
fprintf(1,'\tdim = %d\n', anno.dim);
fprintf(1,'\t#points = %d\n', anno.npts);
fprintf(1,'\tclass = %s\n', anno.ccls);
fprintf(1,'\tkd_ptr = %lu\n', anno.kd_ptr);

author: CHLNDDEV


text:
function  [CQcc, LogP_absCQT, TimeVec, FreqVec, Ures_LogP_absCQT, Ures_FreqVec] = ...
    cqcc(x, fs, B, fmax, fmin, d, cf, ZsdD)

%   Constant Q cepstral coefficients
%   Usage:  CQcc = cqcc(x, fs, B, fmax, fmin, d, cf, ZsdD)
%
%   Input parameters:
%         x        : input signal
%         fs       : sampling frequency
%         B        : number of bins per octave [default = 96]
%         fmax     : highest frequency to be analyzed [default = Nyquist frequency]
%         fmin     : lowest frequency to be analyzed [default = ~20Hz to fullfill an integer number of octave]
%         d        : number of uniform samples in the first octave [default 16]
%         cf       : number of cepstral coefficients excluding 0'th coefficient [default 19]
%         ZsdD     : any sensible combination of the following  [default ZsdD]:
%                      'Z'  include 0'th order cepstral coefficient
%                      's'  include static coefficients (c)
%                      'd'  include delta coefficients (dc/dt)
%                      'D'  include delta-delta coefficients (d^2c/dt^2)
%
%   Output parameters:
%         CQcc              : constant Q cepstral coefficients (nCoeff x nFea)
%         LogP_absCQT       : log power magnitude spectrum of constant Q trasform
%         TimeVec           : time at the centre of each frame [sec]
%         FreqVec           : center frequencies of analysis filters [Hz]
%         Ures_LogP_absCQT  : uniform resampling of LogP_absCQT
%         Ures_FreqVec      : uniform resampling of FreqVec [Hz]
%
%   See also:  cqt
%
%
%   References:
%     M. Todisco, H. Delgado, and N. Evans. A New Feature for Automatic
%     Speaker Verification Anti-Spoofing: Constant Q Cepstral Coefficients.
%     Proceedings of ODYSSEY - The Speaker and Language Recognition
%     Workshop, 2016.
%
%     C. Sch�rkhuber, A. Klapuri, N. Holighaus, and M. D�fler. A Matlab
%     Toolbox for Efficient Perfect Reconstruction log-f Time-Frequecy
%     Transforms. Proceedings AES 53rd Conference on Semantic Audio, London,
%     UK, Jan. 2014. http://www.cs.tut.fi/sgn/arg/CQT/
%
%     G. A. Velasco, N. Holighaus, M. D�fler, and T. Grill. Constructing an
%     invertible constant-Q transform with non-stationary Gabor frames.
%     Proceedings of DAFX11, Paris, 2011.
%
%     N. Holighaus, M. D�fler, G. Velasco, and T. Grill. A framework for
%     invertible, real-time constant-q transforms. Audio, Speech, and
%     Language Processing, IEEE Transactions on, 21(4):775-785, April 2013.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (C) 2016 EURECOM, France.
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International
% License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/4.0/
% or send a letter to
% Creative Commons, 444 Castro Street, Suite 900,
% Mountain View, California, 94041, USA.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Authors: Massimiliano Todisco {todisco [at] eurecom [dot] fr}
%          Hector Delgado {delgado [at] eurecom [dot] fr}
%
% Version: 1.0
% Date: 22.01.16
%
% User are requested to cite the following paper in papers which report
% results obtained with this software package.
%
%     M. Todisco, H. Delgado, and N. Evans. A New Feature for Automatic
%     Speaker Verification Anti-Spoofing: Constant Q Cepstral Coefficients.
%     Proceedings of ODYSSEY - The Speaker and Language Recognition
%     Workshop, 2016.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% CHECK INPUT PARAMETERS
if nargin < 2
    warning('Not enough input arguments.'), return
end

%%% DEFAULT INPUT PARAMETERS
if nargin < 3; B = 96; end
if nargin < 4; fmax = fs/2; end
if nargin < 5; oct = ceil(log2(fmax/20)); fmin = fmax/2^oct; end
if nargin < 6; d = 16; end
if nargin < 7; cf = 19; end
if nargin < 8; ZsdD = 'ZsdD'; end
gamma = 228.7*(2^(1/B)-2^(-1/B));

%%% CQT COMPUTING
Xcq = cqt(x, B, fs, fmin, fmax, 'rasterize', 'full', 'gamma', gamma);

%%% LOG POWER SPECTRUM
absCQT = abs(Xcq.c);
TimeVec = (1:size(absCQT,2))*Xcq.xlen/size(absCQT,2)/fs;
FreqVec = fmin*(2.^((0:size(absCQT,1)-1)/B));
LogP_absCQT = log(absCQT.^2 + eps);

%%% UNIFORM RESAMPLING
kl = (B*log2(1+1/d));
[Ures_LogP_absCQT, Ures_FreqVec] = resample(LogP_absCQT,...
    FreqVec,1/(fmin*(2^(kl/B)-1)),1,1,'spline');

%%% DCT
CQcepstrum = dct(Ures_LogP_absCQT);

%%% DYNAMIC COEFFICIENTS
if strfind(ZsdD, 'Z'); scoeff = 1; else scoeff = 2; end
CQcepstrum_temp = CQcepstrum(scoeff:cf+1,:);
f_d = 3; % delta window size
if strcmp(strrep(ZsdD,'Z',''), 'sdD')
    CQcc = [CQcepstrum_temp; Deltas(CQcepstrum_temp,f_d); ...
        Deltas(Deltas(CQcepstrum_temp,f_d),f_d)];
elseif strcmp(strrep(ZsdD,'Z',''), 'sd')
    CQcc = [CQcepstrum_temp; Deltas(CQcepstrum_temp,f_d)];
elseif strcmp(strrep(ZsdD,'Z',''), 'sD')
    CQcc = [CQcepstrum_temp; Deltas(Deltas(CQcepstrum_temp,f_d),f_d)];
elseif strcmp(strrep(ZsdD,'Z',''), 's')
    CQcc = CQcepstrum_temp;
elseif strcmp(strrep(ZsdD,'Z',''), 'd')
    CQcc = Deltas(CQcepstrum_temp,f_d);
elseif strcmp(strrep(ZsdD,'Z',''), 'D')
    CQcc = Deltas(Deltas(CQcepstrum_temp,f_d),f_d);
elseif strcmp(strrep(ZsdD,'Z',''), 'dD')
    CQcc = [Deltas(CQcepstrum_temp,f_d); Deltas(Deltas(CQcepstrum_temp,f_d),f_d)];
end
end

function D = Deltas(x,hlen)

% Delta and acceleration coefficients
%
% Reference:
%   Young S.J., Evermann G., Gales M.J.F., Kershaw D., Liu X., Moore G., Odell J., Ollason D.,
%   Povey D., Valtchev V. and Woodland P., The HTK Book (for HTK Version 3.4) December 2006.

win = hlen:-1:-hlen;
xx = [repmat(x(:,1),1,hlen),x,repmat(x(:,end),1,hlen)];
D = filter(win, 1, xx, [], 2);
% D = D(:,hlen+1:(end - hlen));
D = D(:,hlen*2+1:end);
D = D./(2*sum((1:hlen).^2));
end

author: nesl


text:
function kel = cel2kel(cel)
%CEL2KEL Convert temperature from celsius to kelvin
%
%  kel = CEL2KEL(cel) converts temperature from celsius to kelvin.
%
%  See also KEL2CEL, CEL2FAR, CEL2RANK.

kel = cel + 273.15;


text:
%MTRAJ Multi-axis trajectory between two points
%
% [Q,QD,QDD] = MTRAJ(TFUNC, Q0, QF, M) is a multi-axis trajectory (MxN) varying
% from configuration Q0 (1xN) to QF (1xN) according to the scalar trajectory function
% TFUNC in M steps. Joint velocity and acceleration can be optionally returned as
% QD (MxN) and QDD (MxN) respectively.  The trajectory outputs have one row per
% time step, and one column per axis.
%
% The shape of the trajectory is given by the scalar trajectory function TFUNC
%      [S,SD,SDD] = TFUNC(S0, SF, M);
% and possible values of TFUNC include @lspb for a trapezoidal trajectory, or
% @tpoly for a polynomial trajectory.
%
% [Q,QD,QDD] = MTRAJ(TFUNC, Q0, QF, T) as above but T (Mx1) is a time
% vector which dictates the number of points on the trajectory.
%
% Notes::
% - If no output arguments are specified Q, QD, and QDD are plotted.
% - When TFUNC is @tpoly the result is functionally equivalent to JTRAJ except
%   that no initial velocities can be specified. JTRAJ is computationally a little
%   more efficient.
%
% See also JTRAJ, MSTRAJ, LSPB, TPOLY.

% Copyright (C) 1993-2015, 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

function [S,Sd,Sdd] = mtraj(tfunc, q0, qf, M)

    if ~isa(tfunc, 'function_handle')
        error('first argument must be a function handle');
    end

    M0 = M;
    if ~isscalar(M)
        M = length(M);
    end
    if numcols(q0) ~= numcols(qf)
        error('must be same number of columns in q0 and qf')
    end

    s = zeros(M, numcols(q0));
    sd = zeros(M, numcols(q0));
    sdd = zeros(M, numcols(q0));

    for i=1:numcols(q0)
        % for each axis
        [s(:,i),sd(:,i),sdd(:,i)] = tfunc(q0(i), qf(i), M);
    end

% - If no output arguments are specified S, SD, and SDD are plotted
%   against time.

    switch nargout
        case 0
            clf

            if isscalar(M0)
                t = [1:M0]';
            else
                t = M0;
            end
            subplot(311)
            plot(t, s); grid; ylabel('s');

            subplot(312)
            plot(t, sd); grid; ylabel('sd');

            subplot(313)
            plot(t, sdd); grid; ylabel('sdd');
            if ~isscalar(M0)
                xlabel('time')
            else
                for c=get(gcf, 'Children');
                    set(c, 'XLim', [1 M0]);
                end
            end
            shg
        case 1
            S = s;
        case 2
            S = s;
            Sd = sd;
        case 3
            S = s;
            Sd = sd;
            Sdd = sdd;
    end

author: Allopart

repo: rbpf-gmapping


text:
function centroids = kMeansInitCentroids(X, K)
%KMEANSINITCENTROIDS This function initializes K centroids that are to be
%used in K-Means on the dataset X
%   centroids = KMEANSINITCENTROIDS(X, K) returns K initial centroids to be
%   used with the K-Means on the dataset X
%

% You should return this values correctly
centroids = zeros(K, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should set centroids to randomly chosen examples from
%               the dataset X
%

randidx = randperm(size(X, 1));
centroids = X(randidx(1:K), :);

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

end

author: zsiciarz

repo: ml-coursera


text:
function max_index_last = i4vec_max_index_last ( n, a )

%*****************************************************************************80
%
%% I4VEC_MAX_INDEX returns the index of the last largest entry in an I4VEC.
%
%  Discussion:
%
%    An I4VEC is a vector of integer values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer A(N), the vector to be searched.
%
%    Output, integer MAX_INDEX_LAST, the index of the last largest entry.
%
  if ( n <= 0 )

    max_index_last = 0;

  else

    amax = a(1);
    max_index_last = 1;

    for i = 2 : n

      if ( amax <= a(i) )
        amax = a(i);
        max_index_last = i;
      end

    end

  end

  return
end