text:
%
% FILE NAME:    miscModelCurationScript_20181012.m
%
% PURPOSE: Script to restore two reactions that were previously deleted,
%          but after further investigation were decided that they should be
%          returned to the model. These reactions were deleted as part of
%          the 'miscModelCurationScript_20180921.m' script, so all of their
%          information was retrieved from the model prior to applying that
%          script, and is re-inserted into the model with this script.
%
%          The reactions are added back to the model in such a way as to
%          try and maintain their original indexing in the model, just for
%          convenience.
%

%% Load model

% load latest version of humanGEM
load('humanGEM.mat');  % v0.4.1

%% Restore transport reactions that were removed earlier

% There was insufficient evidence to remove the reaction that allowed
% transport of HMG-CoA through the peroxisomal membrane (HMGCOAtx).
% Furthermore, the HMG-CoA transport through the mitochondrial membrane
% should proceed via the carnitine shuttle, but for now the simplified
% direct transport is a suitable placeholder reaction (HMR_1572).
%   HMR_1572:  HMG-CoA[c] <=> HMG-CoA[m]
%   HMGCOAtx:  HMG-CoA[c] <=> HMG-CoA[p]

% try to restore them to their original index location in the model (before
% reactions HMR_1284 and HMGCOARr)
[~,r_ind] = ismember({'HMR_1284';'HMGCOARr'},ihuman.rxns);
nMet = length(ihuman.mets);
nGene = length(ihuman.genes);
nProt = length(ihuman.proteins);
met_ind = getIndexes(ihuman,{'HMG-CoA[c]';'HMG-CoA[m]';'HMG-CoA[p]'},'metscomps');

% generate columns for stoich matrix
s1 = zeros(nMet,1);
s1(met_ind(1:2)) = [-1;1];
s2 = zeros(nMet,1);
s2(met_ind([1,3])) = [-1;1];

% restore information in all model fields
ihuman.rxns = [ihuman.rxns(1:r_ind(1)); {'HMR_1572'}; ihuman.rxns((r_ind(1)+1):r_ind(2)); {'HMGCOAtx'}; ihuman.rxns((r_ind(2)+1):end)];
ihuman.S = [ihuman.S(:,1:r_ind(1)), s1, ihuman.S(:,(r_ind(1)+1):r_ind(2)), s2, ihuman.S(:,(r_ind(2)+1):end)];
ihuman.lb = [ihuman.lb(1:r_ind(1)); -1000; ihuman.lb((r_ind(1)+1):r_ind(2)); -1000; ihuman.lb((r_ind(2)+1):end)];
ihuman.ub = [ihuman.ub(1:r_ind(1)); 1000; ihuman.ub((r_ind(1)+1):r_ind(2)); 1000; ihuman.ub((r_ind(2)+1):end)];
ihuman.rev = double((ihuman.lb < 0) & (ihuman.ub > 0));
ihuman.c = [ihuman.c(1:r_ind(1)); 0; ihuman.c((r_ind(1)+1):r_ind(2)); 0; ihuman.c((r_ind(2)+1):end)];
ihuman.rxnNames = [ihuman.rxnNames(1:r_ind(1)); {''}; ihuman.rxnNames((r_ind(1)+1):r_ind(2)); {'Hydroxymethylglutaryl Coenzyme A Reversible Peroxisomal Transport'}; ihuman.rxnNames((r_ind(2)+1):end)];
ihuman.rxnComps = [ihuman.rxnComps(1:r_ind(1)); 3; ihuman.rxnComps((r_ind(1)+1):r_ind(2)); 1; ihuman.rxnComps((r_ind(2)+1):end)];
ihuman.grRules = [ihuman.grRules(1:r_ind(1)); {''}; ihuman.grRules((r_ind(1)+1):r_ind(2)); {''}; ihuman.grRules((r_ind(2)+1):end)];
ihuman.rxnGeneMat = [ihuman.rxnGeneMat(1:r_ind(1),:); zeros(1,nGene); ihuman.rxnGeneMat((r_ind(1)+1):r_ind(2),:); zeros(1,nGene); ihuman.rxnGeneMat((r_ind(2)+1):end,:)];
ihuman.subSystems = [ihuman.subSystems(1:r_ind(1)); {'Transport, mitochondrial'}; ihuman.subSystems((r_ind(1)+1):r_ind(2)); {'Transport, peroxisomal'}; ihuman.subSystems((r_ind(2)+1):end)];
ihuman.eccodes = [ihuman.eccodes(1:r_ind(1)); {''}; ihuman.eccodes((r_ind(1)+1):r_ind(2)); {''}; ihuman.eccodes((r_ind(2)+1):end)];
ihuman.rxnKEGGID = [ihuman.rxnKEGGID(1:r_ind(1)); {''}; ihuman.rxnKEGGID((r_ind(1)+1):r_ind(2)); {''}; ihuman.rxnKEGGID((r_ind(2)+1):end)];
ihuman.rxnEHMNID = [ihuman.rxnEHMNID(1:r_ind(1)); {''}; ihuman.rxnEHMNID((r_ind(1)+1):r_ind(2)); {''}; ihuman.rxnEHMNID((r_ind(2)+1):end)];
ihuman.rxnBiGGID = [ihuman.rxnBiGGID(1:r_ind(1)); {'HMGCOAtm'}; ihuman.rxnBiGGID((r_ind(1)+1):r_ind(2)); {''}; ihuman.rxnBiGGID((r_ind(2)+1):end)];
ihuman.rxnHepatoNET1ID = [ihuman.rxnHepatoNET1ID(1:r_ind(1)); {''}; ihuman.rxnHepatoNET1ID((r_ind(1)+1):r_ind(2)); {''}; ihuman.rxnHepatoNET1ID((r_ind(2)+1):end)];
ihuman.rxnREACTOMEID = [ihuman.rxnREACTOMEID(1:r_ind(1)); {''}; ihuman.rxnREACTOMEID((r_ind(1)+1):r_ind(2)); {''}; ihuman.rxnREACTOMEID((r_ind(2)+1):end)];
ihuman.rxnReferences = [ihuman.rxnReferences(1:r_ind(1)); {''}; ihuman.rxnReferences((r_ind(1)+1):r_ind(2)); {'PMID:14713247'}; ihuman.rxnReferences((r_ind(2)+1):end)];
ihuman.rxnFrom = [ihuman.rxnFrom(1:r_ind(1)); {'HMRdatabase'}; ihuman.rxnFrom((r_ind(1)+1):r_ind(2)); {'Recon3D'}; ihuman.rxnFrom((r_ind(2)+1):end)];
ihuman.rxnMiriams = [ihuman.rxnMiriams(1:r_ind(1)); {''}; ihuman.rxnMiriams((r_ind(1)+1):r_ind(2)); struct('name',{{'pmid'}},'value',{{'PMID:14713247'}}); ihuman.rxnMiriams((r_ind(2)+1):end)];
ihuman.rxnConfidenceScores = [ihuman.rxnConfidenceScores(1:r_ind(1)); NaN; ihuman.rxnConfidenceScores((r_ind(1)+1):r_ind(2)); 0; ihuman.rxnConfidenceScores((r_ind(2)+1):end)];
ihuman.rxnRecon3DID = [ihuman.rxnRecon3DID(1:r_ind(1)); {'HMGCOAtm'}; ihuman.rxnRecon3DID((r_ind(1)+1):r_ind(2)); {''}; ihuman.rxnRecon3DID((r_ind(2)+1):end)];
ihuman.prRules = [ihuman.prRules(1:r_ind(1)); {''}; ihuman.prRules((r_ind(1)+1):r_ind(2)); {''}; ihuman.prRules((r_ind(2)+1):end)];
ihuman.rxnProtMat = [ihuman.rxnProtMat(1:r_ind(1),:); zeros(1,nProt); ihuman.rxnProtMat((r_ind(1)+1):r_ind(2),:); zeros(1,nProt); ihuman.rxnProtMat((r_ind(2)+1):end,:)];
ihuman.priorCombiningGrRules = [ihuman.priorCombiningGrRules(1:r_ind(1)); {''}; ihuman.priorCombiningGrRules((r_ind(1)+1):r_ind(2)); {''}; ihuman.priorCombiningGrRules((r_ind(2)+1):end)];

% clear intermediate variables
clear met_ind nGene nMet nProt r_ind s1 s2

%% Remove non-standard fields (e.g. rxnKEGGID) to comply with defined RAVEN structure
% The following five fields are non-standard. To assist convenient model
% manipulation, they are removed here given that these information have
% been intactly stored elsewhere.
f = {'rxnKEGGID';'rxnEHMNID';'rxnBiGGID';'rxnHepatoNET1ID';'rxnREACTOMEID'};
ihuman = rmfield(ihuman, f);

%% Save model
save('../../model/Human-GEM.mat','ihuman');

author: SysBioChalmers

repo: Human-GEM


text:
%This Matlab script can be used to reproduce Figure 1.12 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Define the range of BS antennas
Mvalues = [1 10 100];

%Angle of the desired UE
varphiDesired = pi/6;

%Range of angles of the interfering UE
varphiInterfererDegrees = -180:1:180;
varphiInterfererRadians = varphiInterfererDegrees*(pi/180);

%Define the antenna spacing (in number of wavelengths)
antennaSpacing = 1/2; %Half wavelength distance

%Preallocate matrix for storing the simulation results
gfunction = zeros(length(varphiInterfererDegrees),length(Mvalues));

%% Go through all number of antennas
for m = 1:length(Mvalues)

    %Generate channel response for the desired UE using (1.23)
    hdesired = exp(1i*2*pi*antennaSpacing*sin(varphiDesired)*(0:Mvalues(m)-1)');

    %Go through all angles of interfering UE
    for n = 1:length(varphiInterfererRadians)

        %Generate channel response for the interfering UE using (1.23)
        hinterfering = exp(1i*2*pi*antennaSpacing*sin(varphiInterfererRadians(n))*(0:Mvalues(m)-1)');

        %Compute the g-function in (1.28), using its definition
        gfunction(n,m) = abs(hdesired'*hinterfering)^2/Mvalues(m);

    end

end

%% Plot the simulation results
figure;
hold on; box on;

plot(varphiInterfererDegrees,gfunction(:,1),'k-','LineWidth',1);
plot(varphiInterfererDegrees,gfunction(:,2),'r--','LineWidth',1);
plot(varphiInterfererDegrees,gfunction(:,3),'b-.','LineWidth',1);

xlabel('Angle of interfering UE [degree]');
ylabel('$$g(\varphi_0^0,\varphi_1^0)$$','Interpreter','latex');

set(gca,'Yscale','log');
xlim([-180 180]);
ylim([1e-5 1e2]);

legend('M=1','M=10','M=100','Location','NorthWest');

author: emilbjornson

repo: massivemimobook


text:
function [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = cf_sum_to_ss(cf2ss)
%% CF_SUM_TO_SS - Sum of several state space models
%
% Syntax:
%   [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = cf_sum_to_ss(cf2ss)
%
% In:
%   cf2ss       - Cell vector of function handles (see below)
%
% Out:
%   F           - Feedback matrix
%   L           - Noise effect matrix
%   Qc          - Spectral density of white noise process w(t)
%   H           - Observation model matrix
%   Pinf        - Covariance of the stationary process
%   dF          - Derivatives of F w.r.t. parameters
%   dQc         - Derivatives of Qc w.r.t. parameters
%   dPinf       - Derivatives of Pinf w.r.t. parameters
%   params      - Input and output parameter information
%
% Description:
%   Calculate the state space formulation corresponding the the following
%   covariance function:
%
%     k(t) = k_1(t) + k_2(t) + ... + k_n(t),
%
%   that is, the state space model of a sum of several state space
%   models. The state space model corresponding to covaraince function
%   k_j is given by the function handle in cf2ss{j} such that:
%
%     [F,L,Qc,H,Pinf,dF,dQc,dPinf] = cf2ss{k}();
%
% Copyright:
%   2012-2014   Arno Solin
%   2013-2014   Jukka Koskenranta
%
%   This software is distributed under the GNU General Public
%   License (version 3 or later); please refer to the file
%   License.txt, included with the software, for details.
%
% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

%% Sums of two state space models by stacking the models

  % Initialize model matrices
  F      = [];
  L      = [];
  Qc     = [];
  H      = [];
  Pinf   = [];
  dF     = [];
  dQc    = [];
  dPinf  = [];
  params = {};

  % Initialize
  params.stationary = true;

  % Stack state space models
  for k=1:length(cf2ss)

    % Form ss model j
    if nargout > 5
        [jF,jL,jQc,jH,jPinf,jdF,jdQc,jdPinf,jparams] = cf2ss{k}();
    else
        [jF,jL,jQc,jH,jPinf] = cf2ss{k}();
    end

    % Stack matrices
    F  = blkdiag(F,jF);
    L  = blkdiag(L,jL);
    Qc = blkdiag(Qc,jQc);
    H  = [H jH];
    Pinf = blkdiag(Pinf,jPinf);

    % Stack derivative matrices
    if nargout > 5

      dF  = mblk(dF,jdF);
      dQc = mblk(dQc,jdQc);
      dPinf = mblk(dPinf,jdPinf);

      % Parameters
      params.stationary = params.stationary & jparams.stationary;

    end
  end

end

function C = mblk(A,B)
% Three-dimensional version of the blk function

  % Get sizes
  sA=size(A); sB=size(B);

  % Check if A or B is empty
  Ae = ~any(sA==0); Be = ~any(sB==0);

  % Numel of sizes to 3
  sA = [sA Ae]; sA = sA(1:3);
  sB = [sB Be]; sB = sB(1:3);

  % Assign space for C
  C = zeros(sA+sB);

  % Set values of A if there is any
  if Ae
      C(1:sA(1), 1:sA(2), 1:sA(3)) = A;
  end

  % Set values of B if there is any
  if Be
      C(sA(1)+(1:sB(1)), ...
        sA(2)+(1:sB(2)), ...
        sA(3)+(1:sB(3))) = B;
  end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function i = iqr(x)
%
% Calculate interquartile range without the stats package
%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

xS = sort(x);         % sort along dimensions
N  = size(x,1);
i  = xS(ceil(3*N/4),:) - xS(ceil(N/4),:);

author: ShapeNet

repo: RenderForCNN


text:
function S=makeCoaster(len,rad,nFrame,doSymmetric)

deltaT=len/1000; deltath=pi/1000; th=-pi/2-deltath; j=1; A=zeros(3,100);
while th>-3*pi/2
  A(:,j)=[-len/2+rad*cos(th) rad*sin(th) 1]'; th=th-deltath; j=j+1;
end
for T=-len/2:deltaT:len/2; A(:,j)=[T rad 1]'; j=j+1; end
th=pi/2-deltath;
while th>-pi/2
  A(:,j)=[len/2+rad*cos(th) rad*sin(th) 1]'; th=th-deltath; j=j+1;
end
for T=len/2:-deltaT:-len/2; A(:,j)=[T -rad 1]'; j=j+1; end
A=A(:,1:j-1);

% Bend the bike chain
nPoint=size(A,2);
if doSymmetric
  center=0;
else
  center=len/3;
end
A(3,:)=6*(0.5-abs((A(1,:)-center)/(len/2+rad)));

% Make it evolve over time
speed=50; nSamp=100; first=1;
nPointSamp=length(1:nSamp:nPoint);

for i=1:nFrame
  inter=first:nSamp:first+nPointSamp*nSamp;
  inter=mod(inter+speed-1,nPoint)+1;
  B=A(:,inter);
  first=first+speed;
  if i==1; S=repmat(B,[1,1,nFrame]); continue; end
  S(:,:,i)=B;
  if rand()<0.1; speed=round(50+rand()*50); end
end

author: vrabaud


text:
function y = vl_nncaffepool(x, poolSize, varargin)
%VL_NNCAFFEPOOL(X) - Apply Caffe-style padded pooling
%    Y = VL_NNCAFFEPOOL(X, POOLSIZE) This layer implements Caffe-style pooling,
%    which applies the pooling kernel to the padding, as well as to the input
%    (unlike MatConvNet which does not include the padding itself in the
%    kernel computation). To mimic this effect, we "pre-pad" the input using
%    an empty convolution, then pass to the standard MatConvNet pooling layer.
%    See the docs fo VL_NNPOOL to understand the full pooling operation.
%
%   VL_NNCAFFEPOOL(..., 'option', value, ...) takes the following option:
%
%   `Method`:: 'max'
%     Specify method of pooling. It can be either 'max' (retain max value
%     over the pooling region per channel) or 'avg' (compute the average
%     value over the pooling region per channel).
%
%   `Pad`:: 0
%     The amount of input padding. Input images are padded with zeros
%     by this number of pixels on all sides before the convolution is
%     computed. It can also be a vector [TOP BOTTOM LEFT RIGHT] to
%     specify a different amount of padding in each direction. The
%     size of the pooling filter has to exceed the padding. Unlike standard
%     VL_NNPOOL, this padding is added before the pooling kernel is applied.
%
% Copyright (C) 2018 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  opts.pad = 0 ;
  opts.stride = 1 ;
  opts.method = 'avg' ;
  opts.extraArgs = {'cuDNN'} ;
  [opts, dzdy] = vl_argparsepos(opts, varargin) ;

  padded = vl_nnconv(x, [], [], 'pad', opts.pad, opts.extraArgs{:}) ;

  if isempty(dzdy)
    y = vl_nnpool(padded, poolSize, 'pad', 0, 'stride', opts.stride, ...
                  'method', opts.method, opts.extraArgs{:}) ;
  else
    dzdy1 = vl_nnpool(padded, poolSize, dzdy{1}, 'pad', 0, 'stride', opts.stride, ...
                  'method', opts.method, opts.extraArgs{:}) ;
    dzdy = vl_nnconv(x, [], [], dzdy1, 'pad', opts.pad, opts.extraArgs{:}) ;
    y = dzdy ;
  end

repo: MFT


text:
function [ tickLabels ] = parse_labels( dcm, labels, idx )
% Generate labels for axis ticks.
%
% This is a protected method of the tapas_Huge class. It cannot be called
% from outside the class.
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2020 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

L = size(dcm.c, 2);
R = dcm.n;

% linear connections
a = cell(size(dcm.a));
for rSource = 1:R
    for rTarget = 1:R
        a{rTarget, rSource} = [labels.regions{rSource} ' => ' ...
            labels.regions{rTarget}];
    end
end

% driving connections
c = cell(size(dcm.c));
for lSource = 1:L
    for rTarget = 1:R
        c{rTarget, lSource} = [labels.inputs{lSource} ' => ' ...
            labels.regions{rTarget}];
    end
end

% modulatory connections
b = cell(size(dcm.b));
for l = 1:L
    for rSource = 1:R
        for rTarget = 1:R
            b{rTarget, rSource, l} = [labels.inputs{lSource} ' MOD ' ...
                labels.regions{rSource} ' => ' labels.regions{rTarget}];
        end
    end
end

% nonlinear connections
d = cell(size(dcm.d));
for rMod = 1:size(d,3)
    for rSource = 1:R
        for rTarget = 1:R
            d{rTarget, rSource, rMod} = [labels.regions{rMod} ' MOD ' ...
                labels.regions{rSource} ' => ' labels.regions{rTarget}];
        end
    end
end

% hemodynamic parameters
hemo = cell(R, 2);
for r = 1:R
    hemo{r, 1} = ['transit ' labels.regions{r}];
    hemo{r, 2} = ['decay ' labels.regions{r}];
end

tmp = [a(:); b(:); c(:); d(:); hemo(:); {'epsilon'}];
tickLabels = tmp([idx.clustering(:); idx.homogenous(:)]);

end

author: translationalneuromodeling

repo: tapas


text:
function load_scope(files,refresh_rate,sample_frequency)
close all
shot_size=round(sample_frequency/refresh_rate);
ring=zeros(shot_size*20,1)';

fig = figure('Position',[500,500,800,600],...
        'NumberTitle','off',...
        'Name','Scope',...
        'doublebuffer','on',...
        'HandleVisibility','on',...
	'Renderer', 'openGL');

plot_ref=plot(zeros(10,1));

for i=1:length(files)
    %files{i};
    fid=fopen(files{i},'r');
    fseek(fid,0,1);
    e_o_f=ftell(fid);
    fseek(fid,0,-1);

    while (ftell(fid)<e_o_f)
        tic
        [data,c]=fread(fid,shot_size,'short');
        data=(data/10000)';
        %data=rand(1,length(data));
        ring=[ring data];
        ring(1:length(data))=[];
        start_display(ring,plot_ref);
        stop_time=toc;
        while stop_time<(1/refresh_rate)
            stop_time=toc;

        end
    end
    fclose(fid)

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [Xh, Xl] = rnd_smp_patch(img_path, type, patch_size, num_patch, upscale)

img_dir = dir(fullfile(img_path, type));

Xh = [];
Xl = [];

img_num = length(img_dir);
nper_img = zeros(1, img_num);

for ii = 1:length(img_dir),
    im = imread(fullfile(img_path, img_dir(ii).name));
    nper_img(ii) = prod(size(im));
end

nper_img = floor(nper_img*num_patch/sum(nper_img));

for ii = 1:img_num,
    patch_num = nper_img(ii);
    im = imread(fullfile(img_path, img_dir(ii).name));
    [H, L] = sample_patches(im, patch_size, patch_num, upscale);
    Xh = [Xh, H];
    Xl = [Xl, L];
end

patch_path = ['Training/rnd_patches_' num2str(patch_size) '_' num2str(num_patch) '_s' num2str(upscale) '.mat'];
save(patch_path, 'Xh', 'Xl');

author: thomas-koehler

repo: SupER


text:
%DEMO_MONOTONIC2  Demonstration of monotonicity information in GP Poisson regression
%
%  Description
%    This demonstration shows the benefit of adding monotonicity
%    information in GP Poisson regression of mortality data.
%
%  References
%
%    Riihimäki and Vehtari (2010). Gaussian processes with
%    monotonicity information.  Journal of Machine Learning Research:
%    Workshop and Conference Proceedings, 9:645-652.
%
%    Broffitt, J. D. (1988). Increasing and increasing convex Bayesian
%    graduation. Transactions of the Society of Actuaries, 40(1), 115-48.
%
%  See also
%    DEMO_*, GP_MONOTONIC, DEMO_MONOTONIC1
%
% Copyright (c) 2014 Ville Tolvanen, Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Note that as number of the observations is quite small, it would
% better to sample hyperparameters with MCMC, but for faster demo CCD
% integration with gp_ia is used below

% Mortality rate data from Broffitt (1988)
% N is the number of people insured under a certain policy
% y is the number of insured who died
S = which('demo_monotonic2');
L = strrep(S,'demo_monotonic2.m','demodata/broffit.txt');
d=readtable(L,'Delimiter',';');
x=d.age;
y=d.y;
z=d.N.*(sum(d.y)./sum(d.N));
zt=repmat(1000,size(y)).*(sum(d.y)./sum(d.N));
[xn,nd.xmean,nd.xstd]=normdata(x);

% Basic GP
pl=prior_t();
pm=prior_sqrtt();
cfl=gpcf_linear('coeffSigma2',10^2,'coeffSigma2_prior',[]);
cfs=gpcf_sexp('lengthScale',.4,'lengthScale_prior',pl,'magnSigma2_prior', pm);
lik=lik_poisson();
% Alternative GP models a) sexp b) lin + sexp
gpa = gp_set('lik', lik, 'cf', {cfs}, 'jitterSigma2', 1e-6, 'latent_method', 'EP');
gpb = gp_set('lik', lik, 'cf', {cfl cfs}, 'jitterSigma2', 1e-6, 'latent_method', 'EP');

% 1) sexp
subplot(2,2,1)
gpiaa=gp_ia(gpa,xn,y,'z',z);
% Do predictions for testing purposes
[Eft,Varft]=gp_pred(gpiaa, xn,y,xn,'z',z);
gp_plot(gpiaa,xn,y,'z',z,'zt',zt,'target','mu','normdata',nd)
hold on
plot(d.age,d.y./d.N*1000,'r*')
axis([34 65 0 23])
xlabel('Age')
ylabel('Deaths/1000')
title('sexp')

% 2) sexp + monotonicity
subplot(2,2,3)
opt=optimset('TolX',1e-1,'TolFun',1e-1,'Display','iter');
gpam=gpa;gpam.xv=xn(2:2:end);
gpam=gp_monotonic(gpam,xn,y,'z',z,'nvd', 1, 'optimize', 'on', ...
                  'opt', opt, 'optimf', @fminlbfgs);
gpiaam=gp_ia(gpam,xn,y,'z',z);
% Do predictions for testing purposes
[Eftm,Varftm]=gp_pred(gpiaam, xn,y,xn,'z',z);
gp_plot(gpiaam,xn,y,'z',z,'zt',zt,'target','mu','normdata',nd)
hold on
plot(d.age,d.y./d.N*1000,'r*')
axis([34 65 0 23])
xlabel('Age')
ylabel('Deaths/1000')
title('sexp + monot.')

% 3) lin + sexp
subplot(2,2,2)
gpiab=gp_ia(gpb,xn,y,'z',z);
gp_plot(gpiab,xn,y,'z',z,'zt',zt,'target','mu','normdata',nd)
hold on
plot(d.age,d.y./d.N*1000,'r*')
axis([34 65 0 23])
xlabel('Age')
ylabel('Deaths/1000')
title('lin + sexp')

% 4) lin + sexp + monotonicity
subplot(2,2,4)
opt=optimset('TolX',1e-1,'TolFun',1e-1,'Display','iter');
gpbm=gpb;gpbm.xv=xn(2:2:end);
gpbm=gp_monotonic(gpb,xn,y,'z',z,'nvd', 1, 'optimize', 'on', ...
                  'opt', opt, 'optimf', @fminlbfgs);
gpiabm=gp_ia(gpbm,xn,y,'z',z);
gp_plot(gpiabm,xn,y,'z',z,'zt',zt,'target','mu','normdata',nd)
hold on
plot(d.age,d.y./d.N*1000,'r*')
axis([34 65 0 23])
xlabel('Age')
ylabel('Deaths/1000')
title('lin + sexp + monot.')


text:
%%*******************************************************************
%% randSDP.m : creates random feasible SDP problems of the form:
%%
%% min { <C,X> : A(X) = b, X psd}
%% where A,C,b are random data.
%%
%% [blk,At,C,b] = randSDP(n,m);
%% n = matrix dimension
%% m = number of linear constraints, i.e. dimension of b.
%%***********************************************************************
%% SDPNAL+
%% Copyright (c) 2014 by
%% Liuqin Yang, Defeng Sun, and Kim-Chuan Toh
%%***********************************************************************

   function  [blk,At,C,b] = randSDP(n,m)

   randn('state',0);
   rand('state',0);

   blk{1,1} = 's'; blk{1,2} = n;
%%
%% generate feasible point X0,Z0
%%
   tmp = randn(n); tmp = tmp*tmp';
   X0 = 0.5*(tmp + tmp');
   tmp = randn(n); tmp = tmp*tmp';
   Z0 = 0.5*(tmp + tmp');
%%
%% set up the matrices Ak,C,b
%%
   b = zeros(m,1);
   y0 = randn(m,1);
   AA  = cell(1,m);
   Aty0 = zeros(n,n);
   for k = 1:m;
      tmp = randn(n,n); Ak = 0.5*(tmp+tmp');
      AA{1,k} = Ak;
      b(k) = sum(sum(Ak.*X0));
      Aty0 = Aty0 + y0(k)*Ak;
   end
   C{1} = Z0 + Aty0;
   At = svec(blk,AA);
%%******************************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function out=Mod2MatMul(matr1,matr2)
%out=Mod2MatMul(matr1,matr2)
%Finds a modulo 2 based matrix product of binary matrices (matr1 & matr2)

%Author: Brhanemedhn Tegegne
%
[r1,c1]=size(matr1);
[r2,c2]=size(matr2);
if c1~=r2
    'Non Matching Matrices'
else
    out=zeros(r1,c2);
    for i=1:r1
        for j=1:c2
            for k=1:c1
                out(i,j)=xor(out(i,j),matr1(i,k)*matr2(k,j));
            end
        end
    end
end
end

author: Sable

repo: mcbench-benchmarks


text:
%ZCROSS Zero-crossing detector
%
% IZ = ZCROSS(IM) is a binary image with pixels set where the corresponding
% pixels in the signed image IM have a zero crossing, a positive pixel
% adjacent to a negative pixel.
%
% Notes::
% - Can be used in association with a Lapalacian of Gaussian image to
%   determine edges.
%
% See also ILOG.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

function iz = zcross(im)

    z = zeros([size(im), 4]);
    K = [1 1 0 ;1 1 0; 0 0 0];
    z(:,:,1) = conv2(im, K, 'same');
    K = [0 1 1 ;0 1 1; 0 0 0];
    z(:,:,2) = conv2(im, K, 'same');
    K = [0 0 0; 1 1 0 ;1 1 0];
    z(:,:,3) = conv2(im, K, 'same');
    K = [0 0 0; 0 1 1 ;0 1 1];
    z(:,:,4) = conv2(im, K, 'same');

    maxval = max(z, [], 3);
    minval = min(z, [], 3);

    iz = double( (maxval > 0) & (minval < 0) );

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = (1 / m) * (X' * X);

[U , S , V] = svd(Sigma);

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

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function lasHeaderInfo = read_las_header(lasFilename)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Reads and displays the header information from a lidar data file in    %
%  the American Society for Photogrammetry & Remote Sensing (ASPRS) LAS   %
%  format, as specified in the following document: "LAS Specification     %
%  Version 1.2" ASPRS, April 29, 2008.  It is backwards-compatible,       %
%  with respect to the LAS version (i.e., header information from lidar   %
%  files in LAS 1.0 and 1.1 formats will display fine).  Requires the     %
%  file "coordinate_sys.mat" to display coordinate system info.           %
%                                                                         %
%  Input:                                                                 %
%     lasFilename: Path and filename of LAS file                          %
%                  Example of calling this function:                      %
%                  headerInfo = read_las_header('C:\my_lidar_file.las');  %
%                                                                         %
%  Output:                                                                %
%     lasHeaderInfo: Up to 40x2 cell arrray containing all the info read  %
%                    in from the LAS file header. If coordinate system    %
%                    info is missing, then lasHeaderInfo will be a 36x2   %
%                    cell array.                                          %
%                                                                         %
%  Chris Parrish                                                          %
%  Created:  5/13/2008                                                    %
%  Modified: 2/26/2009: Minor edits in preparation to post on MATLAB      %
%                       Central. Tested in MATLAB v7.0.1.15 (R14) and     %
%                       Octave v3.0.3.                                    %
%            9/2/2011:  Fixed bug that caused Version Major, Version      %
%                       Minor, and Point Data Format IDs to not be        %
%                       read in correctly. Also now displays vertical     %
%                       datum info and units, if polulated in             %
%                       GeoKeyDirectoryTag.                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Open the LAS file
fid=fopen(lasFilename);
if fid == -1
    error('LAS file could not be opened; check filename and path.')
end

% Read the header info. Numeric data is always read in as a double because
% built-in MATLAB funtions typically expect doubles.
lasHeaderInfo(1,1:2) = [{'File Signature'}; ...
    {sscanf(char(fread(fid,4,'char*1=>char*1')'),'%c')}];
lasHeaderInfo(2,1:2) = [{'File Source ID'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(3,1:2) = [{'Global Encoding'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(4,1:2) = [{'Project ID - GUID data 1'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(5,1:2) = [{'Project ID - GUID data 2'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(6,1:2) = [{'Project ID - GUID data 3'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(7,1:2) = [{'Project ID - GUID data 4'}; ...
    {sscanf(char(fread(fid,8,'char*1=>char*1')'),'%c')}];
lasHeaderInfo(8,1:2) = [{'Version Major'}; ...
    {fread(fid,1,'uint8=>double')}];
lasHeaderInfo(9,1:2) = [{'Version Minor'}; ...
    {fread(fid,1,'uint8=>double')}];
lasHeaderInfo(10,1:2) = [{'System Identifier'}; ...
    {sscanf(char(fread(fid,32,'char*1=>char*1')'),'%c')}];
lasHeaderInfo(11,1:2) = [{'Generating Software'}; ...
    {sscanf(char(fread(fid,32,'char*1=>char*1')'),'%c')}];
lasHeaderInfo(12,1:2) = [{'File Creation Day of Year'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(13,1:2) = [{'File Creation Year'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(14,1:2) = [{'Header Size'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(15,1:2) = [{'Offset to point data'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(16,1:2) = [{'Number of variable length records'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(17,1:2) = [{'Point Data Format ID'}; ...
    {fread(fid,1,'uint8=>double')}];
lasHeaderInfo(18,1:2) = [{'Point Data Record Length'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(19,1:2) = [{'Number of point records'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(20,1:2) = [{'Number of points return 1'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(21,1:2) = [{'Number of points return 2'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(22,1:2) = [{'Number of points return 3'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(23,1:2) = [{'Number of points return 4'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(24,1:2) = [{'Number of points return 5'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(25,1:2) = [{'X scale factor'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(26,1:2) = [{'Y scale factor'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(27,1:2) = [{'Z scale factor'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(28,1:2) = [{'X offset'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(29,1:2) = [{'Y offset'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(30,1:2) = [{'Z offset'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(31,1:2) = [{'Max X'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(32,1:2) = [{'Min X'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(33,1:2) = [{'Max Y'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(34,1:2) = [{'Min Y'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(35,1:2) = [{'Max Z'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(36,1:2) = [{'Min Z'}; ...
    {fread(fid,1,'double=>double')}];

% Read the coordinate system info from the variable length records, if
% available
filePos = lasHeaderInfo{14,2};  % Set file position based on header size
nVarRec = lasHeaderInfo{16,2};  % Number of variable length records
if ~(nVarRec > 0) || ~exist('coordinate_sys.mat')
    disp('Note: Cannot display coordinate system info because')
    disp('coordinate_sys.mat file was not found, or the LAS file does')
    disp('not contain any varaible length records.')
else
    load('-mat', 'coordinate_sys')
    for i = 1:nVarRec
        fseek(fid, filePos, 'bof');
        rsvd = fread(fid,1,'uint16=>double');
        userId = sscanf(char(fread(fid,16,'char*1=>char*1')'),'%c');
        recId = fread(fid,1,'uint16=>double');
        lenRec = fread(fid,1,'uint16=>double');
        desc = sscanf(char(fread(fid,32,'char*1=>char*1')'),'%c');
        filePos = filePos + lenRec + 54;
        % If the current variable length record contains the
        % GeoKeyDirectoryTag, then display the georeferencing info. Note:
        % GeoKeyDirectoryTag = 34735.
        if strcmp(userId,'LASF_Projection') && recId == 34735
            nTags = floor(lenRec/2);
            % Jump to the correct position in the file
            newFilePos = filePos - lenRec;
            fseek(fid, newFilePos, 'bof');
            nCodesFound = 0;     % Number of codes found that match tags
            for j = 1:nTags
                tag = fread(fid,1,'uint16=>double');
                % Read through the CS and units codes and see if the tag
                % matches a known code; if so store it
                index = find(projectedCsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Projected Coordinate System'}; ...
                        projectedCsNames(index)];
                end
                index = find(geographicCsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Geographic Coordinate System'}; ...
                        geographicCsNames(index)];
                end
                index = find(vertCsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Vertical Coordinate System'}; ...
                        vertCsNames(index)];
                end
                index = find(angularUnitsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Angular Units'}; ...
                        angularUnitsNames(index)];
                end
                index = find(linearUnitsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Linear Units'}; ...
                        linearUnitsNames(index)];
                end
            end
        end
    end
end

% Display the header info to the screen in a nice format
disp(['Header information from: ' lasFilename])
for i = 1:length(lasHeaderInfo)
    if isnumeric(lasHeaderInfo{i,2})
        screenOut = [lasHeaderInfo{i,1} ': ' num2str(lasHeaderInfo{i,2})];
    elseif isempty(lasHeaderInfo(i,2))
        screenOut = [lasHeaderInfo{i,1} ': 0'];
    else
        screenOut = [lasHeaderInfo{i,1} ': ' lasHeaderInfo{i,2}];
    end
    disp(screenOut)
end

% Close the LAS file
fclose(fid);


text:
function square_hex_grid_test07 ( )

%*****************************************************************************80
%
%% SQUARE_HEX_GRID_TEST07 tests HEX_GRID_LAYERS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 15;

  nodes_per_layer_test = [ ...
    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 21, 41, 81, 101 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_HEX_GRID_TEST07\n' );
  fprintf ( 1, '  For a hexagonal grid of points in a coordinate box,\n' );
  fprintf ( 1, '  given NODES_PER_LAYER, the number of grid points\n' );
  fprintf ( 1, '  along the first layer,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  HEX_GRID_LAYERS computes LAYERS, the number of\n' );
  fprintf ( 1, '  layers.\n' );

  box(1:2,1:2) = [ 1.0, 2.0; ...
                  4.0, 7.0 ];

  box_print_2d ( box );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   NODES  LAYERS\n' );
  fprintf ( 1, '     PER\n' );
  fprintf ( 1, '   LAYER\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    nodes_per_layer = nodes_per_layer_test ( test );
    layers = hex_grid_layers ( nodes_per_layer, box );
    fprintf ( 1, '  %6d  %6d\n', nodes_per_layer, layers );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rfcn_build()
% rfcn_build()
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

% Compile nms_mex
if ~exist('nms_mex', 'file')
  fprintf('Compiling nms_mex\n');

  mex -O -outdir bin ...
      CXXFLAGS="\$CXXFLAGS -std=c++11"  ...
      -largeArrayDims ...
      functions/nms/nms_mex.cpp ...
      -output nms_mex;
end

if ~exist('nms_gpu_mex', 'file')
   fprintf('Compiling nms_gpu_mex\n');
   addpath(fullfile(pwd, 'functions', 'nms'));
   nvmex('functions/nms/nms_gpu_mex.cu', 'bin');
   delete('nms_gpu_mex.o');
end

repo: R-FCN


text:
function [x_best, psi_best, out] = LLM(mapp, lin_sym_solver, x0, options)
% LLM is a Levenberg-Marquardt algorithm for solving systems of
% nonlinear equations :math:`h(x) = 0`, `x` in :math:`R^m`
% using the nonlinear unconstrained minimization :math:`\textrm{min} \psi(x) = 1/2 ||h(x)||^2`
% s.t. `x` in :math:`R^m`.
%
% USAGE:
%
%    [x_best, psi_best, out] = LLM(mapp, lin_sym_solver, x0, options)
%
% INPUTS:
%    mapp:              function handle provides `h(x)` and gradient `h(x)`
%    lin_sym_solver:    function handle for solving the linear system
%    x0:                initial point
%    options:           structure including the parameteres of scheme
%
%                         * .eta - parameter of the scheme
%                         * .MaxNumIter - maximum number of iterations
%                         * .MaxNumMapEval - maximum number of function evaluations
%                         * .MaxNumGmapEval - maximum number of subgradient evaluations
%                         * .TimeLimit - maximum running time
%                         * .epsilon - accuracy parameter
%                         * .x_opt - optimizer
%                         * .psi_opt - optimum
%                         * .adaptive - update lambda adaptively
%                         * .flag_x_error - 1: saves :math:`x_{error}`, 0: do not saves :math:`x_{error}` (default)
%                         * .flag_psi_error - 1: saves :math:`\psi_{error}`, 0: do not saves :math:`\psi_{error}` (default)
%                         * .flag_time - 1: saves :math:`\psi_{error}`, 0: do not saves :math:`\psi_{error}` (default)
%                         * .Stopping_Crit - stopping criterion
%
%                           1. stop if :math:`||grad|| \leq \epsilon`
%                           2. stop if :math:`||nhxk|| \leq \epsilon`
%                           3. stop if `MaxNumIter` is reached
%                           4. stop if `MaxNumMapEval` is reached
%                           5. stop if `MaxNumGmapEval` is reached
%                           6. stop if `TimeLimit` is reached
%                           7. stop if :math:`||grad|| \leq \textrm{max}(\epsilon, \epsilon^2 * ngradx0)`
%                           8. stop if :math:`||nhxk|| \leq \textrm{max}(\epsilon, \epsilon^2 * nhx0)`
%                           9. stop if (default) :math:`||hxk|| \leq \epsilon` or `MaxNumIter` is reached
%
% OUTPUTS:
%    x_best:            the best approximation of the optimizer
%    psi_best:          the best approximation of the optimum
%    out:               structure including more output information
%
%                         * .T - running time
%                         * .Niter - total number of iterations
%                         * .Nmap - total number of mapping evaluations
%                         * .Ngmap - total number of mapping gradient evaluations
%                         * .merit_func - array including all merit function values
%                         * .x_error - relative error :math:`\textrm{norm}(x_k(:)-x_{opt}(:))/\textrm{norm}(x_{opt})`
%                         * .psi_error - relative error :math:`(\psi_k-\psi_{opt})/(\psi_0-\psi_{opt}))`
%                         * .Status - reason of termination
%
% .. REFERENCE:
% .. 1. M. Ahookhosh, F.J. Aragon Artacho, R.M.T. Fleming, V. Phan, Local convergence of Levenberg-Marquardt methods under Holder metric subregularity, Submitted, (2016)
% .. Author: - Masoud Ahookhosh, System Biochemistry Group, Luxembourg Center for System Biomedicine, University of Luxembourg, Luxembourg
%            - Update: July 2017, M. Ahookhosh

format longG ;

% ================ Error messages for input and output =================
if nargin > 4
    error('The number of input arguments is more than what is needed');
elseif nargin < 4
    error('The number of input arguments is not enough');
end;

if isempty(mapp)
    error('the function handle mapp has to be defined');
elseif ~isa(mapp,'function_handle')
    error('mapp should be a function handle');
end

if isempty(lin_sym_solver)
    error('the function handle lin_sym_solver has to be defined');
elseif ~isa(lin_sym_solver,'function_handle')
    error('lin_sym_solver should be a function handle');
end

if isempty(x0)
    error('The starting point x0 has to be defined');
%elseif ~isa(x0,'numeric')
%    error('x0 should be a numeric vector');
end

% =================== initializing the parameters ======================
% ===== user has requested viewing the default values of "options" =====
[eta,epsilon,MaxNumIter,MaxNumMapEval,MaxNumGmapEval, adaptive, ...
    TimeLimit,flag_x_error,flag_psi_error,flag_time,Stopping_Crit] ...
    = Initialization(options);

if ~isa(eta,'numeric') || (eta <= 0)
    error('eta should be numeric and eta in (0,4*delta)');
end

if isfield(options,'x_opt')
    x_opt=options.x_opt;
elseif flag_x_error==1
    error('x_error requires to x_opt be specified');
end

if flag_x_error == 1
    Nxopt      = sqrt(sum(x_opt(:).^2));
    x_error(1) = sqrt(sum((x0(:)-x_opt(:)).^2))/Nxopt;
end

if flag_psi_error == 1
    psi_error(1) = 1;
end

if flag_time == 1
    Time(1) = 0;
end

xk         = x0;
Xk         = xk;
Niter      = 1;
[hxk,ghxk] = mapp(x0);
Nmap       = 1;
Ngmap      = 1;
grad       = ghxk*hxk;
nhxk       = norm(hxk);
nhx0       = nhxk;
ngradx0    = norm(grad);
I          = eye(length(xk));
merit_func = 0.5*nhxk^2;
StopFlag   = 0;

if adaptive == 1
    lambda = 1;
else
    lambda = options.lambda;
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Main body of LLM.m %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T0 = tic;

% ======================= start of the main loop =======================
while ~StopFlag

    muk = max(lambda*nhxk^eta+(1-lambda)*norm(grad)^eta,1e-8);
    Hk  = ghxk*ghxk'+muk*I;
    dk  = lin_sym_solver(Hk,grad);
    xk  = xk+dk;

    Niter       = Niter+1;
    Xk(:,Niter) = xk;
    [hxk,ghxk]  = mapp(xk);
    Nmap        = Nmap+1;
    Ngmap       = Ngmap+1;
    grad        = ghxk*hxk;
    nhxk        = norm(hxk);

    if adaptive == 1
        l = 0.95^Niter;
        if l >= 1e-2
            lambda = 0.95;
        else
            lambda = max(0.95^Niter,1e-10);
        end
    end

    % ================= Gathering output information ===================
    psik              = 0.5*nhxk^2;
    merit_func(Niter) = psik;
    if flag_time == 1
        Time(Niter+1) = toc(T0);
    end

    if flag_x_error == 1
        Nx_opt = norm(x_opt);
        x_error(Niter+1) = sqrt(sum((xk(:)-x_opt(:)).^2))/Nx_opt;
    end

    if flag_psi_error == 1
        psi_error(Niter+1) = (psik-psi_opt)/(psi0-psi_opt);
    end

    % ================== checking stopping criteria ====================
    T = toc(T0);

    [StopFlag,Status] = StopCriterion(grad,nhxk,Niter,Nmap, ...
    Ngmap,MaxNumIter,MaxNumMapEval,MaxNumGmapEval,T,TimeLimit, ...
    epsilon,nhx0,ngradx0,Stopping_Crit);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Status
x_best         = xk;
psi_best       = psik;
out.Xk         = Xk;
out.T          = T;
out.nhx        = nhxk;
out.merit_func = merit_func';
out.Niter      = Niter;
out.Nmap       = Nmap;
out.Ngmap      = Ngmap;
out.Status     = Status;

if flag_x_error == 1
    out.x_error = x_error;
end
if flag_psi_error == 1
    out.psi_error = psi_error;
end
if flag_time == 1
    out.Time = Time;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of LLM.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: opencobra

repo: cobratoolbox


text:
function K = covGaborard(hyp, x, z, i)

% Gabor covariance function with length scales ell and periods p. The covariance
% function is parameterized as:
%
% k(x,z) = h(x-z), h(t) = exp(-sum(t.^2./(2*ell.^2)))*cos(2*pi*sum(t./p)).
%
% The hyperparameters are:
%
% hyp = [ log(ell_1)
%         log(ell_2)
%          ..
%         log(ell_D)
%         log(p_1)
%         log(p_2)
%          ..
%         log(p_D) ]
%
% For more help on design of covariance functions, try "help covFunctions".
%
% Note that covSM implements a weighted sum of Gabor covariance functions, but
% using an alternative (spectral) parameterization.
%
% Copyright (c) by Hannes Nickisch, 2014-09-26.
%
% See also COVFUNCTIONS.M, COVGABORISO, COVSM.M.

if nargin<2, K = '2*D'; return; end                        % report no of params
if nargin<3, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);                                                % dimensionality
ell = exp(hyp(1:D));                                              % length scale
p = exp(hyp(D+1:2*D));                                                  % period

if dg                                              % compute squared distance d2
  d2 = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
    d2 = sq_dist(diag(1./ell)*x');
  else                                                   % cross covariances Kxz
    d2 = sq_dist(diag(1./ell)*x',diag(1./ell)*z');
  end
end

dp = zeros(size(d2));                                % init sum(t)/p computation
if ~dg
  if xeqz                                                 % symmetric matrix Kxx
    for d=1:D, dp = dp + (x(:,d)*ones(1,size(x,1))-ones(n,1)*x(:,d)')/p(d); end
  else                                                   % cross covariances Kxz
    for d=1:D, dp = dp + (x(:,d)*ones(1,size(z,1))-ones(n,1)*z(:,d)')/p(d); end
  end
end

K = exp(-d2/2);
if nargin<4                                                         % covariance
  K = cos(2*pi*dp).*K;
else                                                               % derivatives
  d = mod(i-1,D)+1;
  if dg, K = K*0; end                                     % handle diagonal case
  if xeqz                                                 % symmetric matrix Kxx
    dd = x(:,d)*ones(1,size(x,1))-ones(n,1)*x(:,d)';
  else                                                   % cross covariances Kxz
    dd = x(:,d)*ones(1,size(z,1))-ones(n,1)*z(:,d)';
  end
  if i<=D                                                         % length scale
    K = (dd/ell(d)).^2.*cos(2*pi*dp).*K;
  elseif i<=2*D                                                         % period
    K = 2*pi*dd/p(d).*sin(2*pi*dp).*K;
  else
    error('Unknown hyperparameter')
  end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%% manually pick neurons
axes(ax_selected);
xrange = get(gca, 'xlim');
yrange = get(gca, 'ylim');
title('click pixels out of the field to end');

while true
    [tmp_x, tmp_y] = ginput(1);
    if tmp_x<xrange(1) || tmp_x>xrange(2) || (tmp_y<yrange(1)) || (tmp_y>yrange(2))
        break;
    else

        % compute distance to all neurons
        if isempty(IDs)
            return;
        end
        % find the closest neuron
        dist = (tmp_x- ctr(IDs,2)).^2 + (tmp_y - ctr(IDs,1)).^2;
        [~, ind] = min(dist);
        IDs(ind) = [];

        %% show spatial shapes of the selected neurons
        img_selected = zeros(size(AA,1), 3);
        col = col0;
        for m=1:3
            img_selected(:, m) = sum(bsxfun(@times, AA(:, IDs), mod(col(:, IDs), 2)), 2);
            col = floor(col/2);
        end
        img_selected = neuron.reshape(img_selected, 2);
        img_selected = img_selected/max(img_selected(:))*(2^16);
        img_selected = uint16(img_selected);

        axes(ax_selected);
        imagesc(img_selected);
        axis equal off tight;
        xlim([xmin-gSiz, xmax+gSiz]);
        ylim([ymin-gSiz, ymax+gSiz]);

        %% show temporal traces of the selected neurons
        tmp_C = neuron.C_raw(IDs, :)';
        tmp_C = bsxfun(@times, tmp_C, 1./max(tmp_C, [], 1));
        axes(ax_trace); cla; hold on;
        for m=1:length(IDs)
            plot(tmp_C(:,m), 'color', color_all(IDs(m),:),  'linewidth', 2);
        end
    end

end

author: zhoupc


text:
function [B1,B2] = box_each_element(V,F)
  % BOX_EACH_ELEMENT
  %
  % [B1,B2] = box_each_element(V,F)
  %
  % Inputs:
  %   V by dim list of mesh vertex positions
  %   F by ss list of simplices indexing rows of V
  % Outputs:
  %   B1 by dim list of min corners of boxes
  %   B2 by dim list of max corners of boxes
  %
  % corners of Axis-aligned boxes containing each element
  [B1,B2] = bounds(reshape(V(F,:)',size(V,2),[],size(F,2)),3);
  B1 = B1';
  B2 = B2';

end

author: alecjacobson

repo: gptoolbox


text:
function A = addcr(c, r)
% Copyright (C2010-2013) Yang XIAO, BJTU, April 16, 2010, E-Mail: yxiao@bjtu.edu.cn.
% This program can produce the Parity Check Matrices of DVBS2 codes, while
% Ref. [1] remindes you that these codes are not good as the LDPC codes IEEE802.16e
% for continuous bits' errors.
% Refenrences:
% [1] X. Huang, Y. Xiao, J. Fan and K. Kim, Girth 4 and low minimum weights' problems of LDPC codes in DVB-S2 and solutions,
% Proceedings of the 5th International Conference on Wireless communications, networking and mobile computing table of contents
% Beijing, China, Pages: 601-604, 2009, ISBN:978-1-4244-3692-7
% [2]  Xiao Y., Kim K, "Alternative good LDPC codes for DVB-S2", 9th International Conference on Signal Processing, 2008. ICSP 2008. Beijing, 26-29 Oct. 2008, Page(s): 1959-1962
% [3]  DVB-S2 standard draft ETSI EN 302 307 V1. 1. 1[S]. 2004, 06.

M = length(c);
N = length(r);
A = zeros(M, N);
for m = 1:M
    A(m, :) = r + c(m);
end


text:
function a = polya_fit_s(data, a, weight)
% POLYA_FIT_S   Maximum-likelihood Dirichlet-multinomial (Polya) precision.
%
% POLYA_FIT_S(data,a) returns the MLE (a) for the matrix DATA,
% subject to a constraint on A/sum(A).
% Each row of DATA is a histogram of counts.
% A is a row vector providing the initial guess for the parameters.
%
% POLYA_FIT_S(data,a,weight) returns the MLE where each histogram is weighted.
% WEIGHT is a column vector of numbers in [0,1] (default all ones).
%
% A is decomposed into S*M, where M is a vector such that sum(M)=1,
% and only S is changed by this function.  In other words, A/sum(A)
% is unchanged by this function.
%
% The algorithm is a generalized Newton iteration, described in
% "Estimating a Dirichlet distribution" by T. Minka.

% Written by Tom Minka

show_progress = 0;

s = sum(a);
m = a/s;

row = (rows(a) == 1);
if row
  sdata = row_sum(data);
else
  [K,N] = size(data);
  sdata = col_sum(data);
end

use_weight = (nargin > 2);

e = [];
% generalized Newton algorithm
for iter = 1:10
  old_s = s;
  if row
    if use_weight
      [g,h,c1,c3] = s_derivatives(a, data, sdata, weight);
    else
      [g,h,c1,c3] = s_derivatives(a, data, sdata);
    end
  else
    if use_weight
      [g,h,c1,c3] = s_derivatives(a, data', sdata', weight');
    else
      [g,h,c1,c3] = s_derivatives(a, data', sdata');
    end
  end
  if g > eps
    r = g + s.*h;
    if r >= 0
      % the maximum is infinity
      s = Inf;
    else
      s = s./(1 + g./h./s);
    end
  end
  if g < -eps & c1 > eps
    s = special_case(s, g, h, c1, c3);
  end
  a = s*m;
  if show_progress
    p = polya_logProb(a, data);
    if use_weight
      p = p.*weight;
    end
    e(iter) = sum(p);
  end
  if show_progress & rem(iter,10) == 0
    plot(e)
    drawnow
  end
  if ~isfinite(s) | abs(s - old_s) < 1e-6
    break
  end
end
if show_progress
  plot(e)
end
%e(iter)

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

function s = special_case(s, g, h, c1, c3)

a1 = h.*s.^2 + c1;
a2 = 2*s.^2.*(h.*s + g);
a3 = s.^3.*(2*g + h.*s);
if abs(2*g + h.*s) < 1e-13
  a3 = c3;
end
b = quad_roots(a1, a2, a3);
a = (g./c1).*((s+b)./b).^2;
% 1/s = 1/s - a
s = 1./(1./s - a);

author: jacobeisenstein

repo: SAGE


text:
function value = r8vec_ascends ( n, x )

%*****************************************************************************80
%
%% R8VEC_ASCENDS determines if an R8VEC is (weakly) ascending.
%
%  Example:
%
%    X = ( -8.1, 1.3, 2.2, 3.4, 7.5, 7.5, 9.8 )
%
%    value = TRUE
%
%    The sequence is not required to be strictly ascending.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the array.
%
%    Input, real X(N), the array to be examined.
%
%    Output, logical VALUE, is TRUE if the entries of X ascend.
%
  for i = 1 : n-1
    if ( x(i+1) < x(i) )
      value = 0;
      return
    end
  end

  value = 1;

  return
end


text:
function display(a)
disp('FAST_MR object:')
disp(struct(a))

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% adds deconvolved traces (sp) to files
% by default it optimizes the neuropil coefficients
% max neuropil coefficients set in ops.maxNeurop

function add_deconvolution(ops, db, iplane0)
ops = build_ops3(db, ops);
ops0 = ops;

% ops.deconvType = 'OASIS' (default) or 'L0'. See
% http://www.biorxiv.org/content/early/2017/06/27/156786 for more info.

% warning: ops0.imageRate now represents the TOTAL frame rate of the recording over all planes. This warning will be disabled in a future version.

for i =  1:length(ops.planesToProcess)
    iplane  = ops0.planesToProcess(i);

    fpath = sprintf('%s/F_%s_%s_plane%d_proc.mat', ops.ResultsSavePath, ...
        ops.mouse_name, ops.date, iplane);
    if exist(fpath, 'file')
        load(fpath);
    else
        fpath = sprintf('%s/F_%s_%s_plane%d.mat', ops.ResultsSavePath, ...
            ops.mouse_name, ops.date, iplane);
        dat = load(fpath);
    end

    if isfield(dat, 'dat')
        dat = dat.dat; % just in case...
    end

    % if z-drift computed, apply correction
    for j = 1:length(dat.Fcell)
        Fcell{j}    = dat.Fcell{j};
        FcellNeu{j} = dat.FcellNeu{j};
    end

    % overwrite fields of ops with those saved to file
    if getOr(ops, 'newops', 0)
        ops = addfields(dat.ops, ops0); % overwrite saved ops with new
    else
        ops = addfields(ops0, dat.ops);% overwrite new ops with saved
    end

    % set up options for deconvolution
    ops.imageRate    = getOr(ops0, {'imageRate'}, 30); % total image rate (over all planes)
    ops.sensorTau    = getOr(ops0, {'sensorTau'}, 2); % approximate timescale in seconds
    ops.sameKernel   = getOr(ops0, {'sameKernel'}, 1); % 1 for same kernel per plane, 0 for individual kernels (not recommended)
    ops.maxNeurop    = getOr(ops0, {'maxNeurop'}, Inf); % maximum neuropil coefficient (default no max)
    ops.recomputeKernel    = getOr(ops0, {'recomputeKernel'}, 0);
    ops.deconvNeuropil = getOr(ops0, {'deconvNeuropil'}, 0);   % whether to deconvolve the neuropil as well

    fprintf('Spike deconvolution, plane %d... \n', iplane)

    % construct Ff and Fneu
    Ff = [];
    Fneu = [];
    for j = 1:numel(Fcell)
        Ff   = cat(1, Ff, Fcell{j}');
        Fneu = cat(1,Fneu, FcellNeu{j}');
    end

    ops.fs                  = ops.imageRate/ops.nplanes;
    ops.estimateNeuropil    = getOr(ops0, 'estimateNeuropil', 1);

    if ops.estimateNeuropil
        [sp, ~, coefs,~, sd, ops, baselines] = wrapperDECONV(ops, Ff, Fneu);
    else
        if isfield(dat.stat, 'neuropilCoefficient')
            coefs = [dat.stat.neuropilCoefficient];
        else
            coefs = .8 * ones(size(Ff,2), 1);
        end
    end

    if ops.deconvNeuropil
        ops.estimateNeuropil = 0;
        spNeu = wrapperDECONV(ops, Fneu);
    end

    stat = dat.stat;
    for j = 1:size(Ff,2)
        stat(j).neuropilCoefficient = coefs(j);
    end

    ops.SpkDcnvAlignToPlane = getOr(ops, 'SpkDcnvAlignToPlane', 0);
    if getOr(ops, 'runningBaseline', 0)
        Ff = Ff - bsxfun(@times, Fneu, [stat.neuropilCoefficient]);
        if ops.SpkDcnvAlignToPlane
            ds = (ops.SpkDcnvAlignToPlane - iplane)/ops.nplanes;
            Ff    = register_F(Ff, ds);
        end
        [sp, ~, ~,~, sd, ops, baselines] = wrapperDECONV(ops, Ff);
    end

    for j = 1:size(Ff,2)
        stat(j).noiseLevel          = sd(j);
        stat(j).baseline           = baselines(j);
    end

    nCum = 0;
    for j = 1:length(dat.Fcell)
       dat.sp{j} = sp(nCum + [1:size(dat.Fcell{j},2)], :)';
       if ops.deconvNeuropil
           dat.spNeu{j} = spNeu(nCum + [1:size(dat.Fcell{j},2)], :)';
       end
       nCum = nCum + size(dat.Fcell{j},2);
    end

    % compute neuropil coefficient and run deconvolution on
    % neuropil-corrected trace
    %     stat = run_deconvolution3(ops, dat);

    dat.stat = stat;
    dat.ops = ops;
    dat.Ff = Ff;

    save(fpath,'-struct', 'dat')
end
%

author: cortex-lab


text:
function ecf_value = ned_to_ecf(ned_value, orp_ecf, position_boolean)
%NED_TO_ECF Convert North-East-Down coordinates to ECF
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////
%
% Convert north, east, and down coordinates to ECF (Earth Centered Fixed).
%
% USAGE:
%   ecf_value = ned_ecf(ned_value, orp_ecf)
%
% INPUTS:
%   ned_value - required : north-east-down coordinates                [m, m, m]
%   orp_ecf - required : ecf coordinates of ORP                       [m, m, m]
%   position_boolean - optional : Is ecf_value a position coordinate? [true, false]
%      Unit normal vectors, velocity vectors, etc. do not require ORP offset.
%
% OUTPUTS:
%   ecf_value - required : ecf coordinates of ned_value               [m, m, m]

ecf_value=(ecf_ned_rot_mat(orp_ecf).')*ned_value(:);
if nargin<3||position_boolean % Assume value is a position coordinate, unless told otherwise
    ecf_value=ecf_value+orp_ecf(:);
end

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

author: ngageoint


text:
function [errorCode,ima_stack,metastructs] = pco_sdk_example_stack(imacount,exposure_time,triggermode)
% set variables and grab images to a Matlab array
%
%   [ima_stack] = pco_sdk_example_stack(imacount,exposure_time,triggermode)
%
% * Input parameters :
%    imacount                number of images to grab
%    exposure_time           camera exposure time (default=10ms)
%    triggermode             camera trigger mode (default=AUTO)
%
% * Output parameters :
%    ima_stack               stack with grabbed images
%    errorCode               ErrorCode returned from pco.camera SDK-functions
%
%grab images from a recording pco.dimax camera
%using functions PCO_AddBufferEx and PCO_WaitforBuffer
%
%function workflow
%open camera
%set variables
%start camera
%grab images
%stop camera
%close camera
%

glvar=struct('do_libunload',0,'do_close',0,'camera_open',0,'out_ptr',[],'err',0);

if(~exist('imacount','var'))
 imacount = 10;
end

if(~exist('exposure_time','var'))
 exposure_time = 5;
end

if(~exist('triggermode','var'))
 triggermode = 0;
end

pco_camera_load_defines();
subfunc=pco_camera_subfunction();
subfunc.fh_lasterr(0);

[errorCode,glvar]=pco_camera_open_close(glvar);
pco_errdisp('pco_camera_setup',errorCode);
disp(['camera_open should be 1 is ',int2str(glvar.camera_open)]);
if(errorCode~=PCO_NOERROR)
 commandwindow;
 return;
end

out_ptr=glvar.out_ptr;

try

subfunc.fh_stop_camera(out_ptr);

cam_desc=libstruct('PCO_Description');
set(cam_desc,'wSize',cam_desc.structsize);
[errorCode,~,cam_desc] = calllib('PCO_CAM_SDK', 'PCO_GetCameraDescription', out_ptr,cam_desc);
pco_errdisp('PCO_GetCameraDescription',errorCode);

if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
 subfunc.fh_enable_timestamp(out_ptr,TIMESTAMP_MODE_BINARYANDASCII);
end

subfunc.fh_set_exposure_times(out_ptr,exposure_time,2,0,2);
subfunc.fh_set_triggermode(out_ptr,triggermode);

%if PCO_ArmCamera does fail no images can be grabbed
errorCode = calllib('PCO_CAM_SDK', 'PCO_ArmCamera', out_ptr);
if(errorCode~=PCO_NOERROR)
 pco_errdisp('PCO_ArmCamera',errorCode);
 ME = MException('PCO_ERROR:ArmCamera','Cannot continue script with not armed camera');
 subfunc.fh_lasterr(errorCode);
 throw(ME);
end

subfunc.fh_set_transferparameter(out_ptr);
subfunc.fh_show_frametime(out_ptr);
ts_bin=subfunc.fh_is_binary_timestamp_enabled(out_ptr);

%get images
[errorCode,ima_stack,metadata_stack,glvar]=pco_camera_stack(imacount,glvar);

subfunc.fh_stop_camera(out_ptr);

if(errorCode==0)
 [~,~,count]=size(ima_stack);
 if(~isempty(metadata_stack))
  libmeta=libstruct('PCO_METADATA_STRUCT');
  [metadatasize,count_m]=size(metadata_stack);
  metastructs(count_m)=get(libmeta);
  for n=1:count_m
   [errorCode_m,metastructs(n)]=subfunc.fh_get_struct_metadata(metadata_stack(:,n),metadatasize);
   if(errorCode_m)
    pco_errdisp('get_struct_metadata',errorCode_m);
    metadata_stack=[];
    break;
   end
  end
  clear libmeta;
 end

 if(count==1)
  m=max(max(ima_stack(10:end-10,10:end-10)));
  disp(['image done maxvalue: ',int2str(m)]);
  if(ts_bin)
   txt=subfunc.fh_print_timestamp_t(ima_stack,1,16);
   disp(['Timestamp data of image: ',txt]);
  end
  if(~isempty(metadata_stack))
   txt=subfunc.fh_print_meta_timestamp(metastructs(1));
   disp(['Timestamp of meta (',num2str(1,'%04d'),'):  ',txt]);
   subfunc.fh_print_meta_struct(metastructs(1));
  end
 else
  disp([int2str(count),' images done']);
  if((ts_bin)||(~isempty(metadata_stack)))
   reply = input('Show timestamps? Y/N [Y]: ', 's');
   if((isempty(reply))||(reply(1)== 'Y')||(reply(1)== 'y'))
    for n=1:count
     if(ts_bin)
      txt=subfunc.fh_print_timestamp_t(ima_stack(:,:,n),1,16);
      disp(['Timestamp data of image(',num2str(n,'%04d'),'): ',txt]);
     end
     if(~isempty(metadata_stack))
      txt=subfunc.fh_print_meta_timestamp(metastructs(n));
      disp(['Timestamp data of meta (',num2str(n,'%04d'),'): ',txt]);
     end
    end
   end
  end
 end
end

catch ME
 errorCode=subfunc.fh_lasterr();
 txt=blanks(101);
 txt=calllib('PCO_CAM_SDK','PCO_GetErrorTextSDK',pco_uint32err(errorCode),txt,100);

 clearvars -except ME glvar errorCode txt;

 if(glvar.camera_open==1)
  glvar.do_close=1;
  glvar.do_libunload=1;
  pco_camera_open_close(glvar);
 end

 if strfind(ME.identifier,'PCO_ERROR:')
  msg=[ME.identifier,' ',ME.message];
  disp(txt);
  warning('off','backtrace')
  warning(msg)
  for k=1:length(ME.stack)
    disp(['from file ',ME.stack(k).file,' at line ',num2str(ME.stack(k).line)]);
  end
  close();
  clearvars -except errorCode;
  commandwindow;
  return;
 else
  close();
  clearvars -except ME;
  rethrow(ME)
 end
end

clearvars -except glvar errorCode ima_stack metastructs;

if(glvar.camera_open==1)
 glvar.do_close=1;
 glvar.do_libunload=1;
 pco_camera_open_close(glvar);
end

if(exist('metastructs','var'))
 [~,count]=size(metastructs);
 if(count>3)
  subfunc=pco_camera_subfunction();
  n=3;
  txt=subfunc.fh_print_meta_timestamp(metastructs(n));
  disp(['Timestamp of meta(',num2str(n,'%04d'),'):  ',txt]);
  subfunc.fh_print_meta_struct(metastructs(n));
  clearvars subfunc;
 end
end

clearvars glvar;
commandwindow;
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function mv = mv_viewMenu(mv,hfig);
%
% mv = mv_viewMenu(mv,hfig);
%
% Add menus with view options for
% the MultiVoxel UI.
%
%
% ras, 04/05
if ieNotDefined('hfig')
    hfig = gcf;
end

if ieNotDefined('mv')
    mv = get(hfig,'UserData');
end

mv.ui.viewMenu = uimenu('ForegroundColor', 'k', 'Label', 'View', 'Separator', 'on');

addFigMenuToggle(mv.ui.viewMenu);

% sort by omniR option
uimenu(mv.ui.viewMenu,'Label','New MultiVoxel UI Window',...
   'Separator','on','Callback','mv_openFig(get(gcf,''UserData'')); multiVoxelUI;');

set(hfig,'UserData',mv);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%PLOTC Plot classifiers
%
%   PLOTC(W,S,LINE_WIDTH)
%   PLOTC(W,LINE_WIDTH,S)
%
% Plots the discriminant as given by the mapping W on predefined axis,
% typically set by scatterd. Discriminants are defined by the points
% where class differences for mapping values are zero.
%
% S is the plot string, e.g. S = 'b--'. In case S = 'col' a color plot is
% produced filling the regions of different classes with different colors.
% Default S = 'k-';
%
% LINE_WIDTH sets the width of the lines and box. Default LINE_WIDTH = 1.5
%
% In W a cell array of classifiers may be given. In S a set of plot strings
% of appropriate size may be given. Automatically a legend is added to
% the plot.
%
% The linear gridsize is read from the global parameter GRIDSIZE, that
% can be set by the function gridsize:  for instance 'gridsize(100)'
% default gridsize is 30. As all these points have to be classified (e.g.
% 100x100) this always done in batch mode. If needed desired the default
% batch size may be changed by PRGLOBAL.
%
% Examples in PREX_CONFMAT, PREX_PLOTC
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, SCATTERD, PLOTM, GRIDSIZE, SETBATCH, PRGLOBAL

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: plotc.m,v 1.6 2009/10/30 11:01:47 davidt Exp $

function handle = plotc(w,varargin)

		linew = 1.5;
	s = [];

	if nargin < 1 | isempty(w)
		handle = prmapping(mfilename,'combiner',s);
		return
	end

	% Extract the plotwidth and linewidth:
	for j = 1:nargin-1
		if isstr(varargin{j})
			s = varargin{j};
		else
			linew = varargin{j};
		end
	end

	ss = char('k-','r-','b-','m-','k--','r--','b--','m--');
	ss = char(ss,'k-.','r-.','b-.','m-.','k:','r:','b:','m:');

	% When we want to plot a list of classifiers, we set up different
	% plot strings ss and call plotc several times.
	if iscell(w)
		w = w(:);
		n = length(w);
		names = [];
		% Check if sufficient plotstrings are available
		if ~isempty(s)
			if size(s,1) == 1
				s = repmat(s,n,1);
			elseif size(s,1) ~= n
				error('Wrong number of plot strings')
			end
		else
			s = ss(1:n,:);
		end
		% Plot the individual boundaries by calling 'mfilename' (i.e. this
		% function again)
		names = [];
		hh = [];
		for i=1:n
			h = feval(mfilename,w{i},deblank(s(i,:)),linew);
			hh = [hh h(1)];
			names = char(names,getname(w{i}));
		end
		% Finally fix the legend:
		names(1,:) = [];
		legend(hh,names,0);
		if nargout > 0
			handle = hh;
		end
		return
	end

	% Now the task is to plot a single boundary (multiple boundaries are
	% already covered above):
	if ~isa(w,'prmapping') | ~istrained(w)
		error('Trained classifier expected')
	end
  %w = w*setbatch;     % Avoid memory prolems with large gridsizes

	[k,c] = size(w);
	c = max(c,2);
	if nargin < 2 | isempty(s)  % default plot string
		s = 1;
	end
	%DXD: Stop when the classifier is not in 2D:
	if (k~=2)
		error('Plotc can only plot classifiers operating in 2D.');
	end

	if ~isstr(s)
		if s > 16 | s < 1
			error('Plotstring undefined')
		else
			s = deblank(ss(s,:));
		end
	end
	% Get the figure size from the current figure:
	hold on
	V=axis;
	hh = [];
	set(gca,'linewidth',linew)

	% linear discriminant
	if isaffine(w) & c == 2 & ~strcmp(s,'col')	% plot as vector
		d = +w;
		n = size(d.rot,2);
		if n == 2, n = 1; end
		for i = 1:n
			w1 = d.rot(:,i); w0 = d.offset(i);
			J = find(w1==0);
			if ~isempty(J)
				w1(J) = repmat(realmin,size(J));
			end
			x = sort([V(1),V(2),(-w1(2)*V(3)-w0)/w1(1),(-w1(2)*V(4)-w0)/w1(1)]);
            if (x(1)==x(2))  % for exactly vertical lines...
                y = V(3:4);
            else
    			y = (-w1(1)*x-w0)/w1(2);
            end
			h = plot(x,y,s);
			set(h,'linewidth',linew)
			hh = [hh h];
		end
	else    % general case: find contour(0)
		% First define the mesh grid:
		n = gridsize;
		m = (n+1)*(n+1);
		dx = (V(2)-V(1))/n;
		dy = (V(4)-V(3))/n;
		[X Y] = meshgrid(V(1):dx:V(2),V(3):dy:V(4));
		D = double([X(:),Y(:),zeros(m,k-2)]*w);
		if min(D(:)) >=0, D = log(D+realmin); end  % avoid infinities

		% A two-class output can be given in one real number, avoid this
		% special case and fix it:
		if c == 2 & min(size(D)) == 1; D = [D -D]; end
		c = size(D,2);

		if ~strcmp(s,'col')

			% Plot the contour lines
			if c < 3
				Z = reshape(D(:,1) - D(:,2),n+1,n+1);
				if ~isempty(contourc([V(1):dx:V(2)],[V(3):dy:V(4)],Z,[0 0]))
					[cc h] = contour([V(1):dx:V(2)],[V(3):dy:V(4)],Z,[0 0],s);
					set(h,'linewidth',linew)
					%DXD Matlab 7 has different handle definitions:
					if str2num(version('-release'))>13,
						h = get(h,'children');
					end
					hh = [hh;h];
				end
			else
				for j=1:c-1
					L = [1:c]; L(j) = [];
					Z = reshape( D(:,j) - max(D(:,L),[],2),n+1,n+1);
					if ~isempty(contourc([V(1):dx:V(2)],[V(3):dy:V(4)],Z,[0 0]))
						[cc h] = contour([V(1):dx:V(2)],[V(3):dy:V(4)],Z,[0 0],s);
						set(h,'linewidth',linew)
						%DXD Matlab 7 has different handle definitions:
						if str2num(version('-release'))>13,
							h = get(h,'children');
						end
						hh = [hh;h];
					end
				end
			end
		else
			% Fill the areas with some colour:
			col = 0; mapp = hsv(c+1); h = [];
			for j=1:c
				L = [1:c]; L(j) = [];
				Z = reshape( D(:,j) - max(D(:,L)',[],1)',n+1,n+1);
				Z = [-inf*ones(1,n+3);[-inf*ones(n+1,1),Z,-inf*ones(n+1,1)];-inf*ones(1,n+3)];
				col = col + 1;
				if ~isempty(contourc([V(1)-dx:dx:V(2)+dx],[V(3)-dy:dy:V(4)+dy],Z,[0 0]))
					[cc h] = contour([V(1)-dx:dx:V(2)+dx],[V(3)-dy:dy:V(4)+dy],Z,[0 0]);
					while ~isempty(cc)
						len = cc(2,1);
						h = [h;fill(cc(1,2:len+1),cc(2,2:len+1),mapp(col,:),'FaceAlpha',0.5)];
						%fill(cc(1,2:len+1),cc(2,2:len+1),mapp(col,:),'FaceAlpha',0.5);
						cc(:,1:len+1) = [];
					end
					hh = [hh;h];
				end
			end
		end
	end
	axis(V);

	% Return the handles if they are requested:
	if nargout > 0, handle = hh; end
	hold off
	return

author: marianux

repo: ecg-kit

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


text:
% This script runs program clear, runmetric
% for nine coordinate systems and puts
% the output in runmetric.tst
if exist('runmetric.tst','file')==2
  delete runmetric.tst
end
diary  runmetric.tst
clear, runmetric(@cylin);
clear, runmetric(@sphr);
clear, runmetric(@cone);
clear, runmetric(@elipcyl);
clear, runmetric(@notort);
clear, runmetric(@oblate);
clear, runmetric(@parab);
clear, runmetric(@toroid);
clear, runmetric(@elipsod);
diary off


text:
function eyemap(SigVector,WindowLength,Xresolution,Yresolution)
% EYEMAP Displays the eye diagram of a signal as an image of intensities
% similarly to a digital oscilloscope
%
% USAGE
%   eyemap(SigVector,WindowLength) displays an image containing the eye
%   diagram of the SigVector as an intensity map with an observation window
%   of length WindowLength.
%
%   SigVector must be a 1xN or Nx1 vector with real terms and WindowLength
%   must be a nonnegative scalar, which is a submultiple of the length of
%   SigVector.
%
%   For a complex signal (e. g. optical signals), EYEMAP will use
%   abs(SigVector).^2 to plot the eye diagram.
%
%   If WindowLength is not a submultiple of the length of SigVector, EYEMAP
%   will truncate SigVector in order to have WindowLength = length(SigVector)/n.
%
%   The generated image will have a default resolution of 256x256.
%
%
%   eyemap(SigVector,WindowLength,Xresolution,Yresolution) displays the eye
%   diagram with a resolution given by Xresolution x Yresolution, with
%   Xresolution and Yresolution as nonnegative integers
%
%
%
% Example for EYEMAP
%
% %Generate random binary sequence
%
% n_simb = 1024; %Number of symbols
% ak=round(rand(1,n_simb)); %create random binary sequence
%
% %Generate binary signal using sequence
% n_pt=128; %Number of samples per symbol
% signal=zeros(1,length(ak)*n_pt); %initialize signal variable
% for k=1:length(ak) %Generate polar signal from digital sequence
%     signal((k-1)*n_pt+1:k*n_pt) = 2*(ak(k)-0.5);
% end;
%
% %Filter signal using a gaussian filter
% signal=real(ifft(fft(signal).*exp(-(1:length(signal)).^2/n_simb^2)));
%
% %Generate random gaussian noise
% noise=randn(size(signal))*0.1;
%
% %Visualize eye diagram of signal and noise
% figure
% eyemap(signal+noise,n_pt*2);
%
%
%
% 2007 Ruben S. Luís

%Verify format of input arguments
switch nargin
    case 0
        %No input arguments, run example
        example;
        return
    case {1,3}
        error('EYEMAP: Not enough input arguments');
    case 2
        %Only 2 input args, use default resolution and contrast
        Xresolution = 256;
        Yresolution = 256;
end;

if nargin>4
    error('EYEMAP: Too many input arguments');
end;
%Remaining option with 5 arguments, we have user defined resolution and
%contrast

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check validity of input arguments %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Check if SigVector is real
if ~isreal(SigVector)
    %SigVector is complex,
    SigVector = abs(SigVector).^2;
end;

%Check size of SigVector
if isempty(SigVector)
    %Empty input signal
    error('EYEMAP: empty input signal');
elseif ~isvector(SigVector)|~isnumeric(SigVector)
    %Input signal not a numeric vector
    error('EYEMAP: a numeric vector input signal is required');
end;

%Check if SigVector is a row or column vector
SigVectorSize = size(SigVector);
if SigVectorSize(1) > SigVectorSize(2)
    %SigVector is a column vector. Rotate to have a row vector
    SigVector = SigVector';
    SigVectorSize = SigVectorSize(end:-1:1);
end;

%Check validity of WindowLength
WindowLength = double(WindowLength);
if ~isscalar(WindowLength)|~isnumeric(WindowLength)|WindowLength<1
    error('EYEMAP: Window Length must be a nonzero integer');
end;

%Check length of SigVector with respect to WindowLength
if rem(SigVectorSize(2),WindowLength)
    %WindowLength is not a sub-multiple of the length of SigVector
    SigVector = SigVector(1:end-rem(SigVectorSize(2),WindowLength));
end;

%Check validity of Xresolution
Xresolution = double(Xresolution);
if ~isscalar(Xresolution)|~isnumeric(Xresolution)|Xresolution<1
    error('EYEMAP: Xresolution must be a nonzero integer');
end;

%Check validity of Yresolution
Yresolution=double(Yresolution);
if ~isscalar(Yresolution)|~isnumeric(Yresolution)|Yresolution<1
    error('EYEMAP: Yresolution must be a nonzero integer');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute Intensity Map %
%%%%%%%%%%%%%%%%%%%%%%%%%

NWindow = length(SigVector)/WindowLength; %Number of observation windows

%Interpolate input signal to the resolution of the X axis
relationX = WindowLength/Xresolution; %relation between windowlength and Xresolution
X = linspace(0,length(SigVector),length(SigVector)/relationX); %new unwindowed X-axis
Y = interp1(SigVector,X,'linear'); %new interpolated signal

%Define grid for quantization of input signal
DY = max(Y)-min(Y); %dynamic range of the signal
Y_grid_ini = linspace(min(Y)-(max(Y)-min(Y))/20,...
    max(Y)+(max(Y)-min(Y))/20,Yresolution); %define grid with +- 5% tolerance
dY = Y_grid_ini(2)-Y_grid_ini(1); %width of the grid

%Quantize the input signal with 16 bit
posY = int16(round((Y-Y_grid_ini(1))/dY)); %Quantized vector of the input signal
Y_grid = int16(0:Yresolution-1); %redefine grid with 16 bit integers

IntensityMap = zeros(Yresolution,Xresolution); %Initialize intensity map

%Intensity map will be calculated by incrementing the values of the cells
%in that map that are crossed by the quantized signal.

counter = 1;
MeanIntensityMap = 0;
for cx = 1:Xresolution %Set a counter for the x-axis of the map
    aux_pos = (0:NWindow-1)*Xresolution+cx; %determine the corresponding position in the signal vector
    auxY = posY(aux_pos); %retain the values of the signal that correspond to that position
    for cy = 1:Yresolution %Set a counter for the y-axis of the map
        aux = sum(auxY==Y_grid(cy)); %find if the signal crosses the cell IntensityMap(cy,cx)
        if aux
            IntensityMap(cy,cx) = aux; %Store the result in the map
            %Compute mean of the intensity map for the contrast
            MeanIntensityMap = MeanIntensityMap+IntensityMap(cy,cx);
            counter = counter + 1;
        end;
    end;
end;
%Normalize Map
MeanIntensityMap=64*MeanIntensityMap/(counter)/max(max(IntensityMap));
IntensityMap=64*IntensityMap/max(max(IntensityMap));

%%%%%%%%%%%%%%%%%%%%%%
% Plot intensity map %
%%%%%%%%%%%%%%%%%%%%%%
%Display image
image((0:Xresolution)/Xresolution,Y_grid_ini,IntensityMap)
set(gca,'Ydir','normal');

%Mean of the map
Ntone=round(2*MeanIntensityMap); %define number of tones in the picture
Ntone=Ntone*(Ntone<=64)+64*(Ntone>64);
Cmap=[linspace(0,1,Ntone) ones(1,64-Ntone-1)]'*ones(1,3); %create inverse colomap
colormap(1-Cmap); %Set color map

return

%%%%%%%%%%%
% EXAMPLE %
%%%%%%%%%%%
function example

%
% Example for EYEMAP
%
%   Generate noisy signal and visualize eye diagram using EYEMAP

%Generate random binary sequence
n_simb = 1024; %Number of symbols
ak=round(rand(1,n_simb)); %create random binary sequence

%Generate binary signal using sequence
n_pt=128; %Number of samples per symbol
signal=zeros(1,length(ak)*n_pt); %initialize signal variable
for k=1:length(ak) %Generate polar signal from digital sequence
    signal((k-1)*n_pt+1:k*n_pt) = 2*(ak(k)-0.5);
end;

%Filter signal using a gaussian filter
signal=real(ifft(fft(signal).*exp(-(1:length(signal)).^2/n_simb^2)));

%Generate random gaussian noise
noise=randn(size(signal))*0.1;

%Visualize eye diagram of signal and noise
figure
eyemap(signal+noise,n_pt*2)
return


text:
function z = mpower( x, y )

%   Disciplined convex/geometric programming information for MPOWER (^):
%      The CVX version of the matrix power function Z=X.^Y supports only
%      the case where X and Y are scalars. In such instances, the rules
%      are identical to those outlined in the help for CVX/POWER.

if length( x ) == 1 && length( y ) == 1,
    z = power( x, y );
else
    error( 'Disciplined convex programming error:\n    Matrix powers not permitted.', 1 ); %#ok
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function options = symbset(varargin)
%SYMBSET  Create or alter the options for Optimization with SymBuilder
%
% options = symbset('param1',value1,'param2',value2,...) creates an
% options structure with the parameters 'param' set to their corresponding
% values in 'value'. Parameters not specified will be set to the SymBuilder
% default.
%
% options = symbset(oldopts,'param1',value1,...) creates a copy of the old
% options 'oldopts' and then fills in (or writes over) the parameters
% specified by 'param' and 'value'.
%
% options = symbset() creates an options structure with all fields set to
% symbset defaults.
%
% symbset() prints a list of all possible fields and their function.

%   Copyright (C) 2013 Jonathan Currie (IPL)

% Print out possible values of properties.
numberargs = nargin;
if (nargin == 0) && (nargout == 0)
    printfields();
    return
end
%Names and Defaults
Names = {'solver','display','cbmode','srckeep','use1stDerivs','use2ndDerivs','preallocate','optiOpts'}';
Defaults = {'auto','off','auto','no','yes','yes','yes',[]}';

%Collect Sizes and lowercase matches
m = size(Names,1);
%Create structure with all names and default values
st = [Names,Defaults]'; options = struct(st{:});

% Check we have char or structure input. If structure, insert arguments
i = 1;
while i <= numberargs
    arg = varargin{i};
    if ischar(arg)
        break;
    end
    if ~isempty(arg)
        if ~isa(arg,'struct')
            error('An argument was supplied that wasn''t a string or struct!');
        end
        for j = 1:m
            if any(strcmp(fieldnames(arg),Names{j,:}))
                val = arg.(Names{j,:});
            else
                val = [];
            end
            if ~isempty(val)
                checkfield(Names{j,:},val);
                options.(Names{j,:}) = val;
            end
        end
    end
    i = i + 1;
end

%Check we have even number of args left
if rem(numberargs-i+1,2) ~= 0
    error('You do not have a value specified for each field!');
end

%Go through each argument pair and assign to correct field
expectval = 0; %first arg is a name
while i <= numberargs
    arg = varargin{i};

    switch(expectval)
        case 0 %field
            if ~ischar(arg)
                error('Expected field name to be a string! - Argument %d',i);
            end
            j = find(strcmp(arg,Names) == 1);
            if isempty(j)  % if no matches
                error('Unrecognised parameter %s',arg);
            elseif(length(j) > 1)
                error('Ambiguous parameter %s',arg);
            end
            expectval = 1; %next arg is a value
        case 1
            checkfield(Names{j,:},lower(arg));
            options.(Names{j,:}) = lower(arg);
            expectval = 0;
    end
    i = i + 1;
end

if expectval %fallen off end somehow
    error('Missing value for %s',arg);
end

%If not using 1st derivatives, turn off 2nd too
if(strcmpi(options.use1stDerivs,'no')), options.use2ndDerivs = 'no'; end

function checkfield(field,value)
%Check a field contains correct data type

if isempty(value)
    return % empty matrix is always valid
end

switch lower(field)
    %String
    case {'use1stderivs','use2ndderivs','preallocate','srckeep'}
        if(ischar(value) && (strcmpi(value,'yes') || strcmpi(value,'no')))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be a string with ''yes'' or ''no''',field);
        end

    case 'solver'
        if(ischar(value) && optiSolver(value))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be a valid solver',field);
        end
    case 'cbmode'
        if(ischar(value) && any(strcmpi(value,{'auto','cppad','ccode','mcode'})))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be ''auto'', ''cppad'', ''ccode'', or ''mcode''',field);
        end
    case 'display'
        if(ischar(value) && any(strcmpi(value,{'off','iter','final'})))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be ''off'', ''iter'' or ''final''',field);
        end
    case {'solveropts','optiopts'}
        if(isstruct(value))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be a structure',field);
        end

    otherwise
        valid = false;
        errmsg = sprintf('Unrecognized parameter name ''%s''.', field);
end

if valid
    return;
else %error
    ME = MException('symbset:FieldError',errmsg);
    throwAsCaller(ME);
end

function printfields()
%Print out fields with defaults
fprintf('        solver: [ Optimization Solver {''auto''} ] \n');
fprintf('       display: [ Solver Display Level {''off''}, ''iter'', ''final'' ]\n');
fprintf('      optiOpts: [ OPTI options (via optiset) - overrides above {[]} ] \n');
fprintf('        cbmode: [ NL Problems: Callback Generation Mode {''auto''}, ''cppad'', ''ccode'', ''mcode'' ]\n');
fprintf('       srckeep: [ NL Problems: Keep C/C++ source generated by SymBuilder {''no''} ]\n');
fprintf('  use1stDerivs: [ NL Problems: Use first derivatives {''yes''} (Turns off use2ndDerivs too if no) ] \n');
fprintf('  use2ndDerivs: [ NL Problems: Use second derivatives {''yes''} ] \n');
fprintf('   preallocate: [ NL Problems: Preallocate return matrices/vectors {''yes''} ] \n');
fprintf('\n');

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% VL_KDTREEBUILD  Build randomized kd-tree
%   FOREST = VL_KDTREEBUILD(X) returns a structure FOREST containing
%   the kd-tree indexing the data X. X is a M x N dimensional matrix
%   of class DOUBLE or SINGLE with one data point per column. Use
%   VL_KDTREEQUERY() to query the index.
%
%   VL_KDETREEBUILD() accepts the following options:
%
%   NumTrees:: [1]
%     Specifies the number of semi-randomized KD-trees to add to the
%     forest. Multiple trees can be used to improve approximate
%     nearest-neighbors searches with VL_KDTREEQUERY().
%
%   ThresholdMethod:: [MEDIAN]
%     Specifies the thresholding method used to split the data.  Two
%     methods are supported: split around the median (MEDIAN) or split
%     around the mean (MEAN).
%
%   Verbose::
%     Increases the verbosity level (can be repeated).
%
%   Distance:: [L2]
%     Specifies the type of norm which the KD-Tree should use (L1 or L2)
%
%   The FOREST structure has the following fields:
%
%   FOREST.TREES::
%     A structure array with one element per tree.
%
%   FOREST.DIMENSION::
%     Dimensionality of the indexed data.
%
%   FOREST.NUMDATA::
%     Number of indexed data points.
%
%   The TREE structure has the following fields:
%
%   TREE.NODES::
%      A structure array representing the nodes of the tree.
%
%   TREE.DATAINDEX::
%      A 1 x NUMDATA vector of class UINT32 representing a permutation
%      of the data.
%
%   Nodes are numbered from 1 to NUMNODES.  The NODES structure array
%   has the following fields:
%
%   NODES.LOWERCHILD and NODES.UPPERCHILD::
%      1 x NUMNODES vectors of class INT32. A positive value is the
%      index of the lower/upper child node. A negative value denotes a
%      leaf and is (after negation) is the first or last element plus
%      one of a range of entries in the permutation TREE.DATAINDEX.
%      Such entries are in turn indexes of the data points that belong
%      to that leaf. Typically there is one point per leaf.
%
%   NODES.SPLITDIMENSION and NODES.SPLITTHRESHOLD::
%      1 x NUMNODES vector of class UINT32 and DOUBLE, respectively,
%      with the index of the splitting dimension and the threshold for
%      each node.
%
%   See also: VL_KDTREEQUERY(), VL_HELP().

% Authors: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function errorflag = check_quadcoefficients_vbmc(gp)
%CHECK_QUADCOEFFICIENTS_VBMC Check that the quadratic coefficients are negative.

% Extract integrated basis functions coefficients
D = size(gp.X,2);
Nb = numel(gp.post(1).intmean.betabar);
betabar = zeros(Nb,numel(gp.post));
for s = 1:numel(gp.post)
    betabar(:,s) = gp.post(s).intmean.betabar;
end
% betabar

if gp.intmeanfun == 3
    errorflag = any(betabar(1+D+(1:D),:) >= 0,2)';
elseif gp.intmeanfun == 4
    tril_mat = tril(true(D),-1);
    tril_vec = tril_mat(:);
    z = zeros(D*D,1);
    errorflag = false;
    for b = 1:size(betabar,2)
        beta_mat = z;
        beta_mat(tril_vec) = betabar(1+2*D+(1:D*(D-1)/2),b);
        beta_mat = reshape(beta_mat,[D,D]);
        beta_mat = beta_mat + beta_mat' + diag(betabar(1+D+(1:D),b));
        try
            [~,dd] = chol(-beta_mat);
        catch
            dd = 1;
        end
        % dd
        errorflag = errorflag | dd;
    end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function y = linspaceNDim(d1, d2, n)
%LINSPACENDIM Linearly spaced multidimensional matrix.
%   LINSPACENDIM(d1, d2) generates a multi-dimensional
%   matrix of 100 linearly equally spaced points between
%   each element of matrices d1 and d2.
%
%   LINSPACENDIM(d1, d2, N) generates N points between
%   each element of matrices X1 and X2.
%
%       Example:
%       d1 = rand(3, 2, 4); d2 = rand(size(d1)); n = 1;
%
%       y = linspaceNDim(d1, d2, n) returns a multidimensional matrix y of
%       size (2, 4, 3, 10)
%
%
%   Class support for inputs X1,X2:
%      float: Multidimensional matrix, vector, double, single
%
%       Steeve AMBROISE --> sambroise@gmail.com
%
%   $ Date: 2009/01/29   21:00:00 GMT $
%   $ revised Date: 2009/02/02   18:00:00 GMT $
%       Bug fixed for singleton dimensions that occur when d1 or d2
%           are empty matrix, scalar or vector.
%
%

if nargin == 2
    n = 100;
end
n  = double(n);
d1 = squeeze(d1); d2 = squeeze(d2);

if ndims(d1)~= ndims(d2) || any(size(d1)~= size(d2))
    error('d1 and d2 must have the same number of dimension and the same size'),
end

NDim = ndims(d1);
%%%%%%%% To know if the two first dimensions are singleton dimensions
if NDim==2 && any(size(d1)==1)
    NDim = NDim-1;
    if all(size(d1)==1)
        NDim = 0;
    end
end

pp      = (0:n-2)./(floor(n)-1);

Sum1    = TensorProduct(d1, ones(1,n-1));
Sum2    = TensorProduct((d2-d1), pp);
y = cat(NDim+1, Sum1  + Sum2, shiftdim(d2, size(d1, 1)==1 ));

%%%%% An old function that I wrote to replace the built in Matlab function:
%%%%% KRON

function Z = TensorProduct(X,Y)
%   Z = TensorProduct(X,Y) returns the REAL Kronecker tensor product of X and Y.
%   The result is a multidimensional array formed by taking all possible products
%   between the elements of X and those of Y.
%
%   If X is m-by-n and Y is p-by-q-by-r, then kron(X,Y)
%   is m-by-p-by-n-by-q-by-r.
%
%   X and Y are multidimensional matrices
%   of any size and number of dimensions
%
%   E.g. if X is of dimensions (4, 5, 3) and Y of dimension (3, 1, 7, 4)
%   TensorProduct(X, Y) returns a multidimensional matrix Z of dimensions:
%   (4, 5, 3, 3, 7, 4)
%
%   $ Date: 2001/11/09 10:20:00 GMT $
%
%       Steeve AMBROISE --> sambroise@gmail.com
%

sX=size(X);sY=size(Y);

ndim1=ndims(X);ndim2=ndims(Y);

indperm=[ndim2+1:ndim1+ndim2,1:ndim2];

% to remove all singleton dimensions
Z=squeeze(repmat(X,[ones(1,ndims(X)),sY]).*...
    permute(repmat(Y,[ones(1,ndims(Y)),sX]),indperm));

author: Sable

repo: mcbench-benchmarks


text:
function gX = dexpKernGradX(kern, x1, x2)

% DEXPKERNGRADX Gradient of the double exponential kernel with respect to a
% point x.
%
% FORMAT
% DESC computes the gradient of the double exponential kernel kernel with
% respect to the input positions.
% ARG kern : kernel structure for which gradients are being computed.
% ARG x1 : locations against which gradients are being computed in the form
% of a design matrix.
% RETURN gX : the returned gradients. The gradients are returned in
% a matrix which is numData x numInputs x numData. Where numData is
% the number of data points and numInputs is the number of input
% dimensions in x1.
%
% FORMAT
% DESC computes the gradident of the double exponential kernel with respect
% to the input positions where both the row positions and column positions
% are provided separately.
% ARG kern : kernel structure for which gradients are being computed.
% ARG x1 : row locations against which gradients are being computed in the
% form of a design matrix.
% ARG x2 : column locations against which gradients are being computed in
% the form of a design matrix.
% RETURN gX : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in x1, numData2 is the number of data
% points in x2 and numInputs is the number of input dimensions in x1
% (currently always one).
%
% SEEALSO dexpKernParamInit, kernGradX, dexpKernDiagGradX
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 3
  x2 = x1;
end
[Nx1, dx1] = size(x1);
[Nx2, dx2] = size(x2);
if (dx1 ~= dx2)
    error('The dimension of x1 and x2 must be the same');
end

% Initialisation of the gradient matrix
gX = zeros(Nx1, dx1, Nx2);

% Computation
K = dexpKernCompute(kern, x1, x2);
for i = size(x1, 1)
    gX(i, :, :) = - kern.decay * transpose(sign(repmat(x1(i, :), Nx2, 1)-x2)) ...
        .* repmat(K(i, :), dx1, 1);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [mask_u, mask_v] = get_pos_err_mask_from_design(design)
%GET_POS_ERR_MASK_FROM_DESIGN Gets mask vectors from the position_errors field
%in the given design. The masks are generated by identifying non-zero elements
%in design.position_errors.
%Input:
%   design - Array design.
%Outputs:
%   mask_u - Mask vector for the position errors along the x-axis.
%   mask_v - Mask vector for the position errors along the y-axis.

if ~isfield(design, 'position_errors') || isempty(design.position_errors)
    mask_u = false(1, design.element_count);
    mask_v = false(1, design.element_count);
    return
end

mask_u = design.position_errors(1,:) ~= 0;
if size(design.position_errors, 1) == 1
    mask_v = false(1, design.element_count);
else
    mask_v = design.position_errors(2,:) ~= 0;
end

end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools