text:
function [ ywinrange ] = bz_MultiLFPPlot( lfp,varargin )
%bz_MultiLFPPlot(lfp). This function plots multiple lfp channels from a
%buzcode lfp structure with the channels appropriately spaced.
%
%USAGE
%   figure
%       bz_MultiLFPPlot(lfp)
%
%INPUT
%   lfp         a buzcode lfp structure. channels are assumed to be ordered
%               from top to bottom
%   (optional)
%   'channels'  subset/ordering of channels to plot (0-index a la neuroscope)
%               default is to take the channels as ordered in lfp.channels
%   'timewin'   only plot a subwindow of time
%   'spikes'    a buzcode spikes struct to display spikes above the LFP
%   'sortmetric'metric by which to sort the cells in the raster (eg FR)
%   'cellgroups'{Ngroups} cell array of logical arrays of group members
%   'plotcells' list of cell UIDs to plot (not implemented. sad.)
%   'axhandle'  axes handle in which to put the plot
%   'scaleLFP'  multiplicative factor to scale the y range of LFP
%   'scalespikes' size of spike points (default:5)
%   'spikeside' 'top' (default) or 'bottom'
%   'LFPlabels' text labels for each LFP channel (default: channel number)
%   'LFPmidpoints' (default: evenly spaced)
%   'lfpcolor'
%   'lfpwidth'
%
%
%DLevenstein 2017
%% parse the inputs!
channelsValidation = @(x) assert(isnumeric(x) || strcmp(x,'all'),...
    'channels must be numeric or "all"');
spikedefault.spindices = [nan nan];

% parse args
p = inputParser;
addParameter(p,'channels','all',@isnumeric)
addParameter(p,'timewin',[0 Inf],@isnumeric)
addParameter(p,'spikes',[]) %should have iscellinfo function
addParameter(p,'sortmetric',[])
addParameter(p,'cellgroups',{})
addParameter(p,'axhandle',gca)
addParameter(p,'scaleLFP',1,@isnumeric)
addParameter(p,'scalespikes',5,@isnumeric)
addParameter(p,'plotcells',nan,@isnumeric)
addParameter(p,'spikeside','top')
addParameter(p,'LFPlabels',[])
addParameter(p,'LFPmidpoints',[])
addParameter(p,'lfpcolor','k')
addParameter(p,'lfpwidth',0.5)
parse(p,varargin{:})
timewin = p.Results.timewin;
channels = p.Results.channels;
spikes = p.Results.spikes;
sortmetric = p.Results.sortmetric;
cellgroups = p.Results.cellgroups;
plotcells = p.Results.plotcells;
ax = p.Results.axhandle;
scaleLFP = p.Results.scaleLFP;
scalespikes = p.Results.scalespikes;
spikeside = p.Results.spikeside;
LFPlabels = p.Results.LFPlabels;
lfpmidpoints = p.Results.LFPmidpoints;
lfpcolor = p.Results.lfpcolor;
lfpwidth = p.Results.lfpwidth;

if isempty(spikes)
    spikes = spikedefault;
else
    %Implement raster sorting - cell sort
    if min(spikes.spindices(:,2))==0
        spikes.spindices(:,2) = spikes.spindices(:,2)+1;
    end

    if isempty(sortmetric)
        sortmetric = 1:max(spikes.spindices(:,2));
    end
    [~,cellsort] =sort(sortmetric);

    %Goups
    if ~isempty(cellgroups)
        for gg = 1:length(cellgroups)
            groupsort{gg} = intersect(cellsort,find(cellgroups{gg}),'stable');
        end
        cellsort = [groupsort{:}];

        %Cell sort should now map from CellID in spikes.spindices to their
        %desired order in the raster, such that cellsort(IDX) is the
        %UID of the IDXth cell in the raster
    end

    %Sort the raster - this fails if every cell doesn't have a group
    [~,sortraster] = sort(cellsort);
    %sortraster(UID) is the position of cell UID in the raster
%     if isempty(sortraster)
%         sortraster = 1:max(spikes.spindices(:,2));
%     end
    %This accounts for cells that have no group
    sortraster(end+1:max(spikes.spindices(:,2)))=nan;
    %

    if ~isnan(plotcells)
        temp = nan(size(sortraster));
        temp(plotcells) = sortraster(plotcells);
        sortraster = temp;
    end
    spikes.spindices(:,2) = sortraster(spikes.spindices(:,2));
    spikes.spindices(isnan(spikes.spindices(:,2)),:) = [];
end

%% Channel and time stuff
%Time Window
windex = lfp.timestamps>=timewin(1) & lfp.timestamps<=timewin(2);
%Channel to data array index mapping
if strcmp(channels,'all')
    chindex = 1:length(lfp.channels);
    channels = lfp.channels;
else
    [~,~,chindex] = intersect(lfp.channels,channels,'stable');
end

winspikes = spikes.spindices(:,1)>=timewin(1) & spikes.spindices(:,1)<=timewin(2);
%% Calculate and implement spacing between channels

%Space based on median absolute deviation over entire recording - robust to outliers.
randtimes = randsample(size(lfp.data,1),1000);
channelrange = 12.*mad(single(lfp.data(randtimes,chindex)),1);
if isempty(lfpmidpoints)
    lfpmidpoints = -cumsum(channelrange);
end
lfp.plotdata = (bsxfun(@(X,Y) X+Y,single(lfp.data(windex,chindex)).*scaleLFP,lfpmidpoints));

switch spikeside
    case 'top'
        spikeplotrange = [0 -lfpmidpoints(1)];
    case 'bottom'
        spikeplotrange = lfpmidpoints(end)+[1.5 0.5].*lfpmidpoints(1);
end
spikes.plotdata = spikes.spindices(winspikes,:);
%spikes.plotdata(:,2) = (spikes.plotdata(:,2)./max(spikes.spindices(:,2))).*(diff(spikeplotrange));
spikes.plotdata(:,2) = bz_NormToRange(spikes.plotdata(:,2),spikeplotrange);
%% Do the plot
ywinrange = fliplr(lfpmidpoints([1 end])+1.*[1 -1].*max(channelrange));
if ~isnan(spikes.spindices)
    switch spikeside
        case 'top'
            ywinrange(2) = ywinrange(2)+max([spikes.plotdata(:,2);0]);
        case 'bottom'
            ywinrange(1) = min(spikes.plotdata(:,2));
    end

end

plot(ax,lfp.timestamps(windex),lfp.plotdata,'color',lfpcolor,'linewidth',lfpwidth)
hold on
plot(ax,spikes.plotdata(:,1),spikes.plotdata(:,2),'k.','markersize',scalespikes)
xlabel('t (s)')
ylabel('LFP Channel')
set(ax,'Ytick',fliplr(lfpmidpoints))

if isempty(LFPlabels)
    set(ax,'yticklabels',fliplr(channels))
else
    set(ax,'yticklabels',fliplr(LFPlabels))
end
ylim(ywinrange)
xlim(timewin)
box off

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: this function mimics the movement of the simulation for target
%           points based on an angular movement
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function test_Rotations()

angVec = 0:0.05:4*pi;

for i=1:length(angVec)

    % Declare how open you want the movements to go
    full_ang = pi/3;
    dbl_ang = 2*full_ang;

    % Find associated angle
    ang1 = mod( angVec(i), dbl_ang );
    if ( (ang1 > full_ang ) && (ang1 <= dbl_ang) )
        ang1 = dbl_ang-ang1;
    end

    ang2 = mod( angVec(i), dbl_ang );
    if ( (ang2 > full_ang ) && (ang2 <= dbl_ang) )
        ang2 = dbl_ang-ang2;
    end

    % Read In Pts!
    [xRef,yRef] = read_File_In('corals.vertex');

    plot(xRef,yRef,'b*'); hold on;
    axis([0 1 0 1]);

    % Store Left/Right Points
    xL_Ref = xRef(1:end/2);     yL_Ref = yRef(1:end/2);
    xR_Ref = xRef(end/2+1:end); yR_Ref = yRef(end/2+1:end);

    % Store Values for Centers of Rotation
    %xL = xRef(1);       yL = yRef(1);
    %xR = xRef(end/2+1); yR = yRef(end/2+1);

    % Rotate Geometry
    %[xR_Ref,yR_Ref] = rotate_Geometry(-ang,xR,yR,xR_Ref,yR_Ref);
    %[xL_Ref,yL_Ref] = rotate_Geometry(ang,xL,yL,xL_Ref,yL_Ref);

    % Store Values for Centers of Rotation
    xL_1 = xRef(1);       yL_1 = yRef(1);           % Left side of Left Pair
    xR_1 = xRef(end/4+1); yR_1 = yRef(end/4+1);     % Right side of Left Pair

    xL_2 = xRef(end/2+1);   yL_2 = yRef(end/2+1);   % Left side of Right Pair
    xR_2 = xRef(3*end/4+1); yR_2 = yRef(3*end/4+1); % Right side of Right Pair

    % -> Rotate Geometry <- %
    % LEFT PAIR %
    [xR_Ref_1,yR_Ref_1] = rotate_Geometry(-ang1,xR_1,yR_1,xRef(end/4+1:end/2),yRef(end/4+1:end/2) );
    [xL_Ref_1,yL_Ref_1] = rotate_Geometry(ang1,xL_1,yL_1,xRef(1:end/4),yRef(1:end/4) );
    % RIGHT PAIR %
    [xR_Ref_2,yR_Ref_2] = rotate_Geometry(-ang2,xR_2,yR_2,xRef(3*end/4+1:end),yRef(3*end/4+1:end) );
    [xL_Ref_2,yL_Ref_2] = rotate_Geometry(ang2,xL_2,yL_2,xRef(end/2+1:3*end/4),yRef(end/2+1:3*end/4) );

    % Store New Geometry
    xRef2 = [xL_Ref_1 xR_Ref_1 xR_Ref_2 xL_Ref_2];
    yRef2 = [yL_Ref_1 yR_Ref_1 yR_Ref_2 yL_Ref_2];

    plot(xRef2,yRef2,'k*'); hold on;
    axis([0 1 0 1]);
    pause(0.01);
    clf;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: rotate geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,y] = rotate_Geometry(ang,xC,yC,xRef,yRef)

len = length(xRef);
x = zeros(1,len); y=x;

xRef = xRef - xC;
yRef = yRef - yC;

for i=1:len
   x(i) =  xRef(i)*cos(ang) - yRef(i)*sin(ang);
   y(i) =  xRef(i)*sin(ang) + yRef(i)*cos(ang);
end

x = x + xC;
y = y + yC;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: reads in info from file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x1,y1] = read_File_In(file_name)

filename = file_name;  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f','CollectOutput',1);

fclose(fileID);        %Close the data file.

mat_info = C{1};   %Stores all read in data

%Store all elements in matrix
mat = mat_info(2:end,1:end);

x1 =  mat(:,1); %store xVals
y1 =  mat(:,2); %store yVals

author: nickabattista


text:
function [dnum_bin, counts_per_bin, sum_per_bin, smallest_per_bin, biggest_per_bin, median_per_bin, std_per_bin, median_time_interval] = bin_irregular(dnum, data, binsize, snum, enum, stepsize)
% BIN_IRREGULAR bin an irregularly-sampled timeseries (like earthquake origin times).
%
%    Usage:
%      [dnum_bin, counts_per_bin, sum_per_bin, smallest_per_bin, biggest_per_bin, median_per_bin, std_per_bin, median_time_interval] = bin_irregular(dnum, data, binsize, snum, enum, [stepsize])
%
%    INPUTS:
%      dnum            - irregular spaced date vector in datenum format
%      data            - data values corresponding to dnum samples
%      binsize         - binsize (in days) to use for output series
%      snum            - start datenum (first centre used for output data)
%      enum            - end datenum (last centre used for output data)
%      stepsize        - (optional) Normally bins do not overlap. But if stepsize is set to a value smaller than binsize, bins will overlap.
%
%    OUTPUTS:
%      dnum_bin        	- regular space date vector (centres of bins)
%      counts_per_bin  	- number of values per bin
%      sum_per_bin     	- sum of all values in each bin
%      smallest_per_bin 	- smallest data value in each bin
%      biggest_per_bin 	- biggest data value in each bin
%      median_per_bin  	- median value in each bin
%      std_per_bin     	- standard deviation of all values in each bin
%      median_time_interval - median time interval between values in each bin
%
%    See also:

% AUTHOR: Glenn Thompson
% $Date$
% $Revision$
l1=length(dnum);
l2=length(data);
dnum_bin = [];
counts_per_bin = [];
sum_per_bin = [];
smallest_per_bin = [];
biggest_per_bin = [];
median_per_bin = [];
std_per_bin = [];
median_time_interval = [];

if (l1==l2)
    dnum_bin = snum+binsize/2 : stepsize : enum-binsize/2; % centres of the bins
    for c=1:length(dnum_bin)
	binstart = dnum_bin(c) - binsize/2; % start of this bin
	binend = dnum_bin(c) + binsize/2;   % end of this bin
	i = find(dnum >= binstart & dnum < binend);
	if length(i)>0
		d = data(i);
		thisdnum = dnum(i);
		counts_per_bin(c) = length(d);
		sum_per_bin(c)=nansum(d);
		median_per_bin(c)=nanmedian(d);
		std_per_bin(c) = std(d);
        	smallest_per_bin(c)=min(d);
        	biggest_per_bin(c)=max(d);
		median_time_interval(c)=median(thisdnum(2:end)-thisdnum(1:end-1));
	else
		counts_per_bin(c) = 0;
		sum_per_bin(c) = 0;
		std_per_bin(c) = 0;
		median_per_bin(c) = NaN;
		smallest_per_bin(c)=NaN;
		biggest_per_bin(c)=NaN;
	end
    end
else
    disp('Could not bin - vector lengths dont match');
end

author: geoscience-community-codes

repo: GISMO


text:
function model = yalmip2gurobimex(interfacedata)
% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
Q       = interfacedata.Q;
c       = interfacedata.c;
K       = interfacedata.K;
integer_variables = interfacedata.integer_variables;
binary_variables = interfacedata.binary_variables;
semicont_variables = interfacedata.semicont_variables;
ub      = interfacedata.ub;
lb      = interfacedata.lb;
x0      = interfacedata.x0;
interfacedata.gettime = 0;
n = length(c);

if ~isempty(ub)
    LB = lb;
    UB = ub;
    LB(binary_variables)  = round(LB(binary_variables));
    UB(binary_variables)  = round(UB(binary_variables));
    LB(integer_variables) = round(LB(integer_variables));
    UB(integer_variables) = round(UB(integer_variables));
else
    LB = [];
    UB = [];
end

if options.showprogress;showprogress('Calling GUROBI',options.showprogress);end

if ~isempty(semicont_variables)
    % Bounds must be placed in LB/UB
    [LB,UB,cand_rows_eq,cand_rows_lp] = findulb(F_struc,K,LB,UB);
    F_struc(K.f+cand_rows_lp,:)=[];
    F_struc(cand_rows_eq,:)=[];
    K.l = K.l-length(cand_rows_lp);
    K.f = K.f-length(cand_rows_eq);

    redundant = find(LB<=0 & UB>=0);
    semicont_variables = setdiff(semicont_variables,redundant);

end

SENSE = 1;     % Minimize
C = full(c);   % Must be full
if size(F_struc,1)>0
    B = full(F_struc(:,1));         % Must be full
    A =-F_struc(:,2:end);
else
    B = [];
    A = [];
end

% Optimized code, make a lot of difference when you make this call 10000
% times in a branch and bound setting...
CTYPE = [char(ones(K.f,1)*61); char(ones(K.l,1)*60)];
VARTYPE = char(ones(length(c),1)*67);
VARTYPE(setdiff(integer_variables,semicont_variables)) = 'I';
VARTYPE(binary_variables)  = 'B';  % Should not happen except from bmibnb
VARTYPE(setdiff(semicont_variables,integer_variables)) = 'S';  % Should not happen except from bmibnb
VARTYPE(intersect(semicont_variables,integer_variables)) = 'N';

% Gurobi assumes semi-continuous variables only can take negative values so
% we negate semi-continuous violating this
NegativeSemiVar = [];
if ~isempty(semicont_variables)
    NegativeSemiVar = find(UB(semicont_variables) < 0);
    if ~isempty(NegativeSemiVar)
        temp = UB(semicont_variables(NegativeSemiVar));
        UB(semicont_variables(NegativeSemiVar)) = -LB(semicont_variables(NegativeSemiVar));
        LB(semicont_variables(NegativeSemiVar)) = -temp;
        A(:,semicont_variables(NegativeSemiVar)) = -A(:,semicont_variables(NegativeSemiVar));
        C(semicont_variables(NegativeSemiVar)) = -C(semicont_variables(NegativeSemiVar));
        if ~isempty(x0)
            x0(NegativeSemiVar) = -NegativeSemiVar;
        end
    end
end

if nnz(Q)>0
    [ii,jj,kk] = find(Q);
    ii = ii-1;
    jj = jj-1;
    comp = computer;
    % According to Wotao's testing
    if isempty(strfind(comp,'64')) | isequal(comp,'GLNXA64') | isequal(comp,'PCWIN64')
        options.gurobi.QP.qrow = int32(ii)';
        options.gurobi.QP.qcol = int32(jj)';
    elseif strfind(comp,'MACI64')
        options.gurobi.QP.qrow = int32(ii)';
        options.gurobi.QP.qcol = int32(jj)';
    else
        options.gurobi.QP.qrow = int64(ii)';
        options.gurobi.QP.qcol = int64(jj)';
    end
    options.gurobi.QP.qval = kk';
end

if ~options.verbose
    options.gurobi.DisplayInterval = 0;
    options.gurobi.Display = 0;
end

if ~isempty(x0)
    options.gurobi.Start = x0(:)';
end

if options.savedebug
    save gurobidebug
end

if ~isempty(K.sos.type)
    options.gurobi.SOS.weights = spalloc(length(c),length(K.sos.type),0);
    for i = 1:length(K.sos.type)
        options.gurobi.SOS.types(i)= int32(str2num(K.sos.type(i)));
        options.gurobi.SOS.weights(K.sos.variables{i},i) = K.sos.weight{i};
    end
end
model.C = C;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [inH2O] = Pa2inH2O(Pa)
% Convert pressure from pascals to inches of 4-degree-C water.
% Chad Greene 2012
inH2O = Pa*0.00401463;

author: Sable

repo: mcbench-benchmarks


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

function Tmat=Tsfer(gamma,beta,r)
%%Transformation in spherical coordinates
%%input: alfa angle along z axis, beta angle along y axis, r trasl along z
%%along z axis
%%output: trasformation matrix
Tmat=RotZ(gamma)*RotY(beta)*Tras(0,0,r);

end


text:
classdef MatrixOperator < spx.dict.Operator

    properties(SetAccess=private)
        % The matrix as an operator
        A
        % The Hermitian of matrix
        AH
    end

    methods
        function self = MatrixOperator(A)
            self.A = A;
            self.AH = A';
        end

        function [mm, nn]  = get_size(self)
            [mm, nn] = size(self.A);
        end

        function result = apply(self, vectors)
            result = self.A * vectors;
        end

        function result = apply_transpose(self, vectors)
            result = transpose(self.A) * vectors;
        end

        function result = apply_ctranspose(self, vectors)
            result = self.AH * vectors;
        end

        function result = norm(self)
            result = norm(self.A);
        end

    end

    methods
        % We override some of the implementations here
        function result = apply_columns(self, vectors, columns)
            result = self.A(:, columns) * vectors;
        end

        function result = columns(self, columns)
            result = self.A(:, columns);
        end

        function result = columns_operator(self, columns)
            b = self.A(:, columns);
            result = spx.dict.MatrixOperator(b);
        end

        function result = double(self)
            result = self.A;
        end

        function result = subsref2(self, s)
            result = subsref(self.A, s);
        end
    end

end

author: indigits

repo: sparse-plex


text:
function pred_value = xval_simple_ols_featureselect(X, Y, pthreshold, holdout_method)
% Check: a very simple leave-one-out cross-validated regression
% pred_value = xval_simple_ols_loo(X, Y, p-value selection for univariate feature selection, holdout_method)
%
% pred_value: cross-validated predictions of outcome data
% X: n x variables matrix of predictors
% Y: n x 1 vector of outcomes
%
% Tor Wager, June 2010
%
% Go to any LASSO output directory and run this:
%
% maskInfo = iimg_read_img(fullfile(pwd, 'mask.img'), 2);
% dat{1} =  iimg_get_data(maskInfo, anticimages);
% pred_value = xval_simple_ols_loo(X, Y)

[N, k] = size(X);
pred_value = zeros(N, 1);

holdout_set = nested_select_holdout_set;

create_figure('test', 1, 2);

fprintf('Fold: %03d', 0);
for wh_fold = 1:length(holdout_set)

    fprintf('\b\b\b%3.0f ', wh_fold);

    % select training data
    Xi = X;
    Yi = Y;

    Yi(holdout_set{wh_fold}) = [];
    Xi(holdout_set{wh_fold}, :) = [];              % leave out the missing observation(s)

    Xtest = X(holdout_set{wh_fold}, :);
    Ytest = Y(holdout_set{wh_fold}, :);          % only used later when we test prediction

    ntrain = size(Xi, 1);
    nholdout = size(Xtest, 1);

    % select features based on univariate correlations
    [r, p, Tstat] = correlation_fast_series(Xi, Yi);
    wh_features = p <= pthreshold;
    nfeatures(wh_fold) = sum(wh_features);
    if nfeatures(wh_fold) == 0
        disp('Warning: no features pass threshold');
        wh_features = p <= prctile(p, 10);
    end

    Xi = Xi(:, wh_features);

    Xi = [ones(ntrain, 1) Xi]; % add intercept

    % Make prediction
    b = pinv(Xi) * Yi;
    pred_value(holdout_set{wh_fold}, 1) = [ones(nholdout, 1) Xtest(:, wh_features)] * b;

    create_figure('test', 1, 2, 1); subplot(1, 2, 1); plot(Xi*b, Yi, 'ko');
    plot(pred_value(holdout_set{wh_fold}, 1), Ytest, 'ro', 'MarkerFaceColor', 'r');
    subplot(1, 2, 2);
   title('Black circles = training set; Red = holdout obs');
    drawnow;

end

cm = colormap(jet(N));
figure; hold on;
for i = 1:N
    plot(pred_value(i), Y(i), 'ko', 'MarkerFaceColor', cm(i, :));
end
xlabel('Predicted outcome (xval)'); ylabel('Outcome');
title('Color = order in data series');

function holdout_set = nested_select_holdout_set
        % purpose: return holdout_set variable

        nobs_s = length(Y);

        switch lower(holdout_method)
            case 'loo'
                holdout_set = cell(1, nobs_s);
                for i = 1:nobs_s, holdout_set{i} = i; end

            case 'l4o_covbalanced'
                disp('Selecting holdout sets: Balancing on covariates and outcome, and also trying to ensure that each obs is selected equally often.');
                holdout_proportion = 4 ./ nobs_s; % prop for leave-4-out
                nfolds = 40;
                if isempty(cov_val)
                    wh_holdout = xval_select_holdout_set(Y, [], nfolds, holdout_proportion, verbose);
                else
                    wh_holdout = xval_select_holdout_set(Y, cov_val, nfolds, holdout_proportion, verbose);
                end
                holdout_set = cell(1, nfolds);
                for k = 1:nfolds
                    holdout_set{k} = wh_holdout(:, k);
                end

            case 'categorical_covs'

                holdout_set = xval_select_holdout_set_categoricalcovs(cov_val);

            case 'balanced4'
                nfolds = 4;
                holdout_set = cell(1, nfolds);
                [ys, indx] = sort(Y);
                for k = 1:nfolds

                    holdout_set{k} = indx(k:nfolds:end);
                    if isempty(holdout_set{k}), error('Holdout set construction error: Dataset too small?'); end

                end

            otherwise error('Unknown holdout method.  See help.');
        end
end

end % main function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function scale = p04_scale ( neqn )

%*****************************************************************************80
%
%% P04_SCALE returns scale factors for problem p04.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = 17.7;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [colatitude, longitude, r] = AMICO_Cart2sphere(x, y, z)
	r = sqrt(x.^2 + y.^2 + z.^2);
	longitude  = atan2(y, x);
	colatitude = atan2(sqrt(x.^2 + y.^2), z);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
clear
% close all
clc;

addpath(genpath('../libs'));
path_to_matconvnet = '/home/skong2/scratch/matconvnet-1.0-beta23_modifiedDagnn';

run(fullfile(path_to_matconvnet, 'matlab', 'vl_setupnn'));
addpath(genpath(fullfile('dependencies', 'matconvnet','examples')));
%% read matconvnet model
load('imdb_toydata_v3_from_mnist.mat');
imdb.path = './toydata_v3';
imdb.path_to_dataset = './toydata_v3';

% set GPU
gpuId = 2; %[1, 2];
gpuDevice(gpuId);
flagSaveFig = true; % {true false} whether to store the result

saveFolder = 'main001_instSeg_v1_absEucMM/';
modelName = 'softmax_net-epoch-31.mat';
%% modify network for testing
netMat = load( fullfile('./exp', saveFolder, modelName) );
netMat = netMat.net;
netMat = dagnn.DagNN.loadobj(netMat);

for loopIdx = 5:-1:1
    rmLayerName = sprintf('loop%d_meanshift_Y_l2norm', loopIdx);
    netMat.vars(netMat.layers(netMat.getLayerIndex(rmLayerName)).outputIndexes).precious = 1;

    % add regression loss
    rmLayerName = sprintf('loop%d_instSeg_reg', loopIdx);
    netMat.removeLayer(rmLayerName); % remove layer

    % add max-margin loss
    rmLayerName = sprintf('loop%d_instSeg_MM', loopIdx);
    netMat.removeLayer(rmLayerName); % remove layer

    % add max-margin loss
    rmLayerName = sprintf('loop%d_meanshift_cosSim', loopIdx);
    netMat.removeLayer(rmLayerName); % remove layer
end

rmLayerName = 'obj_instSeg_MM';
if ~isnan(netMat.getLayerIndex(rmLayerName))
    netMat.removeLayer(rmLayerName); % remove layer
end
rmLayerName = 'obj_instSeg_reg';
if ~isnan(netMat.getLayerIndex(rmLayerName))
    netMat.removeLayer(rmLayerName); % remove layer
end
rmLayerName = 'res7_cosSim';
if ~isnan(netMat.getLayerIndex(rmLayerName))
    netMat.removeLayer(rmLayerName); % remove layer
end
netMat.vars(netMat.layers(netMat.getLayerIndex('res7_l2norm')).outputIndexes).precious = 1;
%% configure training environment
saveFolder = [strrep(saveFolder, '/', '') '_visualization'];
if ~isdir(saveFolder) && flagSaveFig
    mkdir(saveFolder);
end
netMat.move('gpu');
% netMat.mode = 'test';
netMat.mode = 'normal';
testList = find(imdb.set==2);
rawCountsAll = {};

semanticMaskMat = [];
instanceMaskMat = [];
imgMat = [];
predInstanceMaskMat_loop0 = [];
predInstanceMaskMat = {};
for loopIdx = 1:5
    predInstanceMaskMat{loopIdx} = [];
end
st = 0;
for i = 1:5%length(testList)
    curBatch = fullfile(imdb.path_to_dataset, imdb.imgList{testList(i)});
    datamat = load(curBatch);

    datamat.semanticMaskMat = reshape(datamat.semanticMaskMat, [size(datamat.semanticMaskMat,1), size(datamat.semanticMaskMat,2), 1, size(datamat.semanticMaskMat,3)]);
    datamat.instanceMaskMat = reshape(datamat.instanceMaskMat, [size(datamat.instanceMaskMat,1), size(datamat.instanceMaskMat,2), 1, size(datamat.instanceMaskMat,3)]);

    imFeed = bsxfun(@minus, datamat.imgMat, imdb.meta.meanvalue);
    inputs = {'input', gpuArray(imFeed)};
    netMat.eval(inputs) ;
    res7_l2norm = gather(netMat.vars(netMat.layers(netMat.getLayerIndex('res7_l2norm')).outputIndexes).value);
    predInstanceMaskMat_loop0(:,:,:,st+1:st+size(datamat.imgMat,4)) = res7_l2norm;

    for loopIdx = 1:5
        curLayerName = sprintf('loop%d_meanshift_Y_l2norm', loopIdx);
        predInstanceMaskMat{loopIdx}(:,:,:,st+1:st+size(datamat.imgMat,4)) = ...
            gather(netMat.vars(netMat.layers(netMat.getLayerIndex(curLayerName)).outputIndexes).value);
    end

    imgMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.imgMat;
    semanticMaskMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.semanticMaskMat;
    instanceMaskMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.instanceMaskMat;
    %predInstanceMaskMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = res7_l2norm;
    fprintf('\t%d/%d\n', i, length(testList));

    st = st + size(datamat.imgMat,4);
end
%% grouping results to visualize
num = 25; % the total number of example to show
M = 5; % the number of rows in the panel, each row shows N images as below
N = 5; % the number of columns in the panel, each column shows M images as above
% h = 28;
% w = 28;
panelSZ = round(64/2)*2; % the size (height/width) of the small square image
% stepSize = (panelSZ-h)/2;

% image
showcaseSamples = stitchImage2panel(imgMat, panelSZ, M, N, 1, num);

% ground-truth semantic mask
showcaseSemanticMask= semanticMaskMat(:,:,:,1:num);
showcaseSemanticMask = reshape(showcaseSemanticMask, [numel(showcaseSemanticMask(:,:,1,1)), num]);
showcaseSemanticMask = showStitch(showcaseSemanticMask, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);

% ground-truth instance mask
showcaseInstanceMask= instanceMaskMat(:,:,:,1:num);
showcaseInstanceMask = reshape(showcaseInstanceMask, [numel(showcaseInstanceMask(:,:,1,1)), num]);
showcaseInstanceMask = showStitch(showcaseInstanceMask, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);

predInstanceMaskMat_loop0 = stitchImage2panel(predInstanceMaskMat_loop0, panelSZ, M, N, 1, num);
predInstanceMaskMat_loop0 = predInstanceMaskMat_loop0 - min(predInstanceMaskMat_loop0(:));
predInstanceMaskMat_loop0 = predInstanceMaskMat_loop0 ./ max(predInstanceMaskMat_loop0(:));

predInstanceMask = {};
for loopIdx = 1:5
    predInstanceMask{loopIdx} = stitchImage2panel(predInstanceMaskMat{loopIdx}, panelSZ, M, N, 1, num);
    predInstanceMask{loopIdx} = predInstanceMask{loopIdx} - min(predInstanceMask{loopIdx}(:));
    predInstanceMask{loopIdx} = predInstanceMask{loopIdx} ./ max(predInstanceMask{loopIdx}(:));
end
%% save the results
figFolder = './figFolder';
if ~isdir(figFolder)
    mkdir(figFolder);
end
a = strfind(modelName,'epoch');
path_to_save = fullfile(figFolder, [strrep(saveFolder,'/','') '_' modelName(a(1):end-4)]);
if ~isdir(path_to_save)
    mkdir(path_to_save);
end

% summary vidualization
imgFig2 = figure(1);
subWindowH = 2;
subWindowW = 2;
set(imgFig2, 'Position', [100 100 1400 900]) % [1 1 width height]
windowID = 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseSamples); axis off image; title('random samples'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseSemanticMask); axis off image; title('showcase SemanticMask'); caxis([0, 11]); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseInstanceMask); axis off image; title('showcase instanceMask'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMaskMat_loop0); axis off image; title('predInstanceMask'); colorbar; windowID = windowID + 1;
if flagSaveFig
    export_fig( sprintf('%s/fig00_visualization.jpg', path_to_save) );
end

imgFig2 = figure(2);
subWindowH = 2;
subWindowW = 3;
set(imgFig2, 'Position', [100 100 1400 900]) % [1 1 width height]
windowID = 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMaskMat_loop0); axis off image; title('loop0'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMask{1}); axis off image; title('loop1'); caxis([0, 11]); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMask{2}); axis off image; title('loop2'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMask{3}); axis off image; title('loop3'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMask{4}); axis off image; title('loop4'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMask{5}); axis off image; title('loop5'); colorbar; windowID = windowID + 1;
if flagSaveFig
    export_fig( sprintf('%s/fig01_visualization_looping.jpg', path_to_save) );
end

if flagSaveFig
    showcaseSemanticMask = uint8(showcaseSemanticMask * 255/11);
    showcaseInstanceMask = showcaseInstanceMask*255/5;

    imwrite(showcaseSamples, fullfile(path_to_save, sprintf('fig02_showcaseSamples.bmp')));
    imwrite(showcaseSemanticMask, fullfile(path_to_save, sprintf('fig03_showcaseSemanticMask.bmp')));
    imwrite(showcaseInstanceMask, fullfile(path_to_save, sprintf('fig04_showcaseInstanceMask.bmp')));
    imwrite(predInstanceMaskMat_loop0, fullfile(path_to_save, sprintf('fig05_predInstanceMask-loop0.bmp')));
    imwrite(predInstanceMask{1}, fullfile(path_to_save, sprintf('fig06_predInstanceMask-loop1.bmp')));
    imwrite(predInstanceMask{2}, fullfile(path_to_save, sprintf('fig07_predInstanceMask-loop2.bmp')));
    imwrite(predInstanceMask{3}, fullfile(path_to_save, sprintf('fig08_predInstanceMask-loop3.bmp')));
    imwrite(predInstanceMask{4}, fullfile(path_to_save, sprintf('fig09_predInstanceMask-loop4.bmp')));
    imwrite(predInstanceMask{5}, fullfile(path_to_save, sprintf('fig10_predInstanceMask-loop5.bmp')));

    save(sprintf('%s/results.mat', path_to_save), 'imgMat', 'instanceMaskMat', 'semanticMaskMat', 'predInstanceMaskMat');
end

%% leaving blank

author: aimerykong


text:
function cdf = beta_cdf ( x, a, b )

%*****************************************************************************80
%
%% BETA_CDF evaluates the Beta CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0D+00 < A,
%    0.0D+00 < B.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x <= 0.0 )
    cdf = 0.0;
  elseif ( x <= 1.0 )
    cdf = beta_inc ( a, b, x );
  else
    cdf = 1.0;
  end

  return
end


text:
function x = gpml_randn(seed,varargin)

% Generate pseudo-random numbers in a quick and dirty way.
% The function makes sure, we obtain the same random numbers using Octave and
% Matlab for the demo scripts.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch 2010-07-07.

if nargin==2
  sz = varargin{1};
else
  sz = zeros(1,nargin-1);
  for i=1:nargin-1
    sz(i) = varargin{i};
  end
end

n = prod(sz); N = ceil(n/2)*2;

% minimal uniform random number generator for uniform deviates from [0,1]
% by Park and Miller
a = 7^5; m = 2^31-1;
% using Schrage's algorithm
q = fix(m/a); r = mod(m,a); % m = a*q+r
u = zeros(N+1,1); u(1) = fix(seed*2^31);
for i=2:N+1
  % Schrage's algorithm for mod(a*u(i),m)
  u(i) = a*mod(u(i-1),q) - r*fix(u(i-1)/q);
  if u(i)<0, u(i) = u(i)+m; end
end
u = u(2:N+1)/2^31;

% Box-Muller transform: Numerical Recipies, 2nd Edition, $7.2.8
% http://en.wikipedia.org/wiki/Box-Muller_transform
w = sqrt(- 2*log(u(1:N/2)));                             % split into two groups
x = [w.*cos(2*pi*u(N/2+1:N)); w.*sin(2*pi*u(N/2+1:N))];
x = reshape(x(1:n),sz);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [params] = evalValidTest(model, validData, testData, params)
  startTime = clock;
  [validCosts] = evalCost(model, validData, params);
  [testCosts] = evalCost(model, testData, params);

  validCounts = initCosts();
  validCounts = updateCounts(validCounts, validData);
  validCosts = scaleCosts(validCosts, validCounts);

  testCounts = initCosts();
  testCounts = updateCounts(testCounts, testData);
  testCosts = scaleCosts(testCosts, testCounts);

  modelStr = wInfo(model);
  endTime = clock;
  timeElapsed = etime(endTime, startTime);

  params.curTestPerpWord = exp(testCosts.word);
  [params.scaleTrainCosts] = scaleCosts(params.trainCosts, params.trainCounts);
  logStr = sprintf('# eval %.2f, %d, %d, %.2fK, %.2f, train=%.2f, valid=%.2f, test=%.2f,%s, time=%.2fs', params.curTestPerpWord, ...
    params.epoch, params.iter, params.speed, params.lr, ...
    params.scaleTrainCosts.total, validCosts.total, testCosts.total, modelStr, timeElapsed);
  fprintf(2, '%s\n', logStr);
  fprintf(params.logId, '%s\n', logStr);

  if validCosts.total < params.bestCostValid
    params.bestCostValid = validCosts.total;
    params.costTest = testCosts.total;
    params.testPerplexity = params.curTestPerpWord;

    fprintf(2, '  save model test perplexity %.2f to %s\n', params.testPerplexity, params.modelFile);
    fprintf(params.logId, '  save model test perplexity %.2f to %s\n', params.testPerplexity, params.modelFile);
    save(params.modelFile, 'model', 'params');

    if params.saveHDF
      saveHDF5([params.modelFile '.h5'], model, params);
    end
  end
end

author: lmthang


text:
% Formatting_step1_cleanup

clear all;close all;clc

code_dir = 'C:\Users\Xiu\Dropbox\Github\FishExplorer';
addpath(genpath(code_dir));

%% Set Manaully!
% file directories
M_dir = GetFishDirectories();
M_tcutoff = {3500,[],[],3500,3600,4000,1800,[],... % Fish 1-8
    [],5000,[] }; % FIsh 9-11

fpsec = 1.97; % Hz

% also check manual frame-correction section below

%%
poolobj=parpool(8);
%%
range_fish = 15;

for i_fish = range_fish,
    tic
    disp(['i_fish = ', num2str(i_fish)]);

    %% load data
    disp(['load data: fish ' num2str(i_fish)]);
    data_dir = M_dir{i_fish};
    if exist(fullfile(data_dir,'cell_resp_dim_lowcut.mat'), 'file'),
        load(fullfile(data_dir,'cell_resp_dim_lowcut.mat'));
    elseif exist(fullfile(data_dir,'cell_resp_dim.mat'), 'file'),
        load(fullfile(data_dir,'cell_resp_dim.mat'));
    else
        errordlg('find data to load!');
    end

    load(fullfile(data_dir,'cell_info.mat'));
    if exist(fullfile(data_dir,'cell_resp_lowcut.stackf'), 'file'),
        cell_resp_full = read_LSstack_fast_float(fullfile(data_dir,'cell_resp_lowcut.stackf'),cell_resp_dim);
    elseif exist(fullfile(data_dir,'cell_resp.stackf'), 'file'),
        cell_resp_full = read_LSstack_fast_float(fullfile(data_dir,'cell_resp.stackf'),cell_resp_dim);
    else
        errordlg('find data to load!');
    end

    load(fullfile(data_dir,'frame_turn.mat'));
    frame_keys = frame_turn;

    %% load anatomy
    tiffname = fullfile(data_dir,'ave.tif');
    info = imfinfo(tiffname,'tiff');
    nPlanes = length(info);
    s1 = info(1).Height;
    s2 = info(1).Width;
    ave_stack = zeros(s1,s2,nPlanes);
    for i=1:nPlanes,
        ave_stack(:,:,i) = imread(tiffname,i);
    end

    %% fix the left-right flip in the anatomy stack and subsequent cell_info
    anat_stack = fliplr(ave_stack);

    [s1,s2,~] = size(ave_stack);
    for i_cell = 1:length(cell_info),
        % fix '.center'
        cell_info(i_cell).center(2) = s2-cell_info(i_cell).center(2)+1; %#ok<*SAGROW>
    end

    %% reformat coordinates
    numcell_full = cell_resp_dim(1);
    temp = [cell_info(:).center];
    XY = reshape(temp',[],numcell_full)';
    Z = [cell_info.slice]';
    CellXYZ = horzcat(XY,Z);

    %% Manual occasional frame corrections: ONLY RUN ONCE!!!!!!!!!

    if true,
        % add 1 frame at start: ------------- for which fish again? all till F11?
        cell_resp_full = horzcat(cell_resp_full(:,1),cell_resp_full);
        frame_keys = vertcat(frame_keys(1,:),frame_keys);
    end

    % correction of an error in Fish #1
    if i_fish == 1,
        cell_resp_full = horzcat(cell_resp_full(:,1:20),cell_resp_full);
    end

    %% Crop end of experiment (when cell-segmentation drift > ~1um, manually determined)
    if ~isempty(M_tcutoff{i_fish}),
        cell_resp_full = cell_resp_full(:,1:M_tcutoff{i_fish});
        frame_keys = frame_keys(1:M_tcutoff{i_fish},:);
    end

    %% detrend
    CR_dtr = zeros(size(cell_resp_full));
    tmax=size(cell_resp_full,2);
    numcell_full=size(cell_resp_full,1);

    parfor i=1:numcell_full,
        cr = cell_resp_full(i,:);
        crd = 0*cr;
        for j=1:100:tmax,
            if j<=150,
                tlim1 = 1;
                tlim2 = 300;
            elseif j>tmax-150,
                tlim1 = tmax-300;
                tlim2 = tmax;
            else
                tlim1 = j-150;
                tlim2 = j+150;
            end
            crr = cr(tlim1:tlim2);
            crd(max(1,j-50):min(tmax,j+50)) = prctile(crr,15);
        end
%         if mod(i,100)==0,
%             disp(num2str(i));
%         end
        CR_dtr(i,:) = cr-crd;
    end

    TimeSeries = single(CR_dtr);

    %% Place holder for registered coordinates from morphing to ZBrain
    CellXYZ_norm = CellXYZ;

    %% Save files

    % directory to save data formatted for distribution:
    save_masterdir = GetNestedDataDir();
    save_dir = fullfile(save_masterdir,['subject_' num2str(i_fish)]);
    if ~exist(save_dir, 'dir')
        mkdir(save_dir);
    end

    % save time-series
    h5create(fullfile(save_dir,'TimeSeries.h5'),'/TimeSeries',size(TimeSeries),'Datatype','single','ChunkSize',[1000 100]);
    h5write(fullfile(save_dir,'TimeSeries.h5'),'/TimeSeries',TimeSeries);

    % save .mat files
    filename = fullfile(save_dir,'CoreInfo.mat');
    varList = {'CellXYZ','anat_stack','fpsec'};
    save(filename,varList{:});

    filename = fullfile(save_dir,'OptionalInfo.mat');
    varList = {'numcell_full','CellXYZ_norm'};
    save(filename,varList{:});

    filename = fullfile(save_dir,'AdditionalInfo.mat');
    save(filename,'frame_keys');

    %%
    toc
end

%%
% delete(poolobj);

author: xiuyechen

repo: FishExplorer


text:
function outd = figmenus(d)
% Set up sliceomatic's gui menus within structure D

% Main Figure Menu
  set(gcf,'menubar','none');

  % File menu
  d.filemenu = uimenu(gcf,'label','File');
  d.fcopy = uimenu(d.filemenu, 'label', 'Copy figure','callback', 'sliceomatic copy');
  d.fprint  = uimenu(d.filemenu,'label','Print...','callback','sliceomatic print');
  % How do get these props onto the print figure?
  %d.fprints = uimenu(d.filemenu,'label','Print Setup...','callback','printdlg -setup');
  % ---
  d.fexit = uimenu(d.filemenu, 'label', 'Close','callback','closereq',...
                   'separator','on');

  % Controls Menu
  d.defcontrols = uimenu(gcf,'label','Controls', 'callback',@controlmenu);
  d.camtoolbar = uimenu(d.defcontrols,'label','Camera toolbar','callback', 'sliceomatic cameratoolbar');
  d.dcalpha = uimenu(d.defcontrols,'label','Controls Transparency');
  d.dcalpha1= uimenu(d.dcalpha,'label','1','callback','sliceomatic controlalpha 1');
  d.dcalpha8= uimenu(d.dcalpha,'label','.8','callback','sliceomatic controlalpha .8');
  d.dcalpha6= uimenu(d.dcalpha,'label','.6','callback','sliceomatic controlalpha .6');
  d.dcalpha5= uimenu(d.dcalpha,'label','.5','callback','sliceomatic controlalpha .5');
  d.dcalpha4= uimenu(d.dcalpha,'label','.4','callback','sliceomatic controlalpha .4');
  d.dcalpha2= uimenu(d.dcalpha,'label','.2','callback','sliceomatic controlalpha .2');
  d.dcalpha0= uimenu(d.dcalpha,'label','0','callback','sliceomatic controlalpha 0');
  d.dclabels= uimenu(d.defcontrols','label','Tick Labels','callback','sliceomatic controllabels');
  d.dcvis   = uimenu(d.defcontrols','label','Visible','callback','sliceomatic controlvisible');
  %  d.dcslice = uimenu(d.defcontrols,'label','Slice Controls','callback','sliceomatic useslicecontrols');
%  d.dciso   = uimenu(d.defcontrols,'label','Iso Surface Control','callback','sliceomatic useisocontrols','separator','on');

  % Remove this once we have more controls to enable and disable.
%  set(d.defcontrols,'vis','off');

  % Default for new slices menu
  d.defmenu = uimenu(gcf,'label','Object_Defaults', 'callback', @defaultmenu);
  d.dfacet  = uimenu(d.defmenu,'label','Slice Color Faceted','callback','sliceomatic defaultfaceted');
  d.dflat   = uimenu(d.defmenu,'label','Slice Color Flat',   'callback','sliceomatic defaultflat');
  d.dinterp = uimenu(d.defmenu,'label','Slice Color Interp', 'callback','sliceomatic defaultinterp');
  d.dtex    = uimenu(d.defmenu,'label','Slice Color Texture','callback','sliceomatic defaulttexture');
  d.dcnone  = uimenu(d.defmenu,'label','Slice Color None','callback','sliceomatic defaultcolornone');
  d.dtnone  = uimenu(d.defmenu,'label','Slice Transparency None','callback','sliceomatic defaulttransnone','separator','on');
  d.dtflat  = uimenu(d.defmenu,'label','Slice Transparency Flat','callback','sliceomatic defaulttransflat');
  d.dtinterp= uimenu(d.defmenu,'label','Slice Transparency Interp','callback','sliceomatic defaulttransinterp');
  d.dttex   = uimenu(d.defmenu,'label','Slice Transparency Texture','callback','sliceomatic defaulttranstexture');
  d.dlflat  = uimenu(d.defmenu,'label','IsoSurface Lighting Flat','callback','sliceomatic defaultlightflat','separator','on');
  d.dlsmooth= uimenu(d.defmenu,'label','IsoSurface Lighting Smooth','callback','sliceomatic defaultlightsmooth');
  d.dcflat  = uimenu(d.defmenu,'label','Contour Color Flat',   'callback','sliceomatic defaultcontourflat','separator','on');
  d.dcinterp= uimenu(d.defmenu,'label','Contour Color Interp', 'callback','sliceomatic defaultcontourinterp');
  d.dcblack = uimenu(d.defmenu,'label','Contour Color Black',  'callback','sliceomatic defaultcontourblack');
  d.dcwhite = uimenu(d.defmenu,'label','Contour Color White',  'callback','sliceomatic defaultcontourwhite');
  d.dclinew = uimenu(d.defmenu,'label','Contour Line Width');
  d.dcl1    = uimenu(d.dclinew,'label','1','callback','sliceomatic defaultcontourlinewidth 1');
  d.dcl2    = uimenu(d.dclinew,'label','2','callback','sliceomatic defaultcontourlinewidth 2');
  d.dcl3    = uimenu(d.dclinew,'label','3','callback','sliceomatic defaultcontourlinewidth 3');
  d.dcl4    = uimenu(d.dclinew,'label','4','callback','sliceomatic defaultcontourlinewidth 4');
  d.dcl5    = uimenu(d.dclinew,'label','5','callback','sliceomatic defaultcontourlinewidth 5');
  d.dcl6    = uimenu(d.dclinew,'label','6','callback','sliceomatic defaultcontourlinewidth 6');

  d.defcolor='texture';
  d.defalpha='none';
  d.deflight='smooth';
  d.defcontourcolor='black';
  d.defcontourlinewidth=1;

  % Set props for all slices menu
  d.allmenu = uimenu(gcf,'label','AllSlices');
  uimenu(d.allmenu,'label','Color Faceted','callback','sliceomatic allfacet');
  uimenu(d.allmenu,'label','Color Flat','callback','sliceomatic allflat');
  uimenu(d.allmenu,'label','Color Interp','callback','sliceomatic allinterp');
  uimenu(d.allmenu,'label','Color Texture','callback','sliceomatic alltex');
  uimenu(d.allmenu,'label','Color None','callback','sliceomatic allnone');
  uimenu(d.allmenu,'label','Transparency None','callback','sliceomatic alltnone','separator','on');
  uimenu(d.allmenu,'label','Transparency .5','callback','sliceomatic alltp5');
  uimenu(d.allmenu,'label','Transparency Flat','callback','sliceomatic alltflat');
  uimenu(d.allmenu,'label','Transparency Interp','callback','sliceomatic alltinterp');
  uimenu(d.allmenu,'label','Transparency Texture','callback','sliceomatic allttex');

  % Set Help menu
  d.Help = uimenu(gcf,'label','Help');
  uimenu(d.Help,'label','About dicomrt_explore','callback','sliceomatic showabout_dicomrt_explore');
  uimenu(d.Help,'label','About sliceomatic','callback','sliceomatic showabout_sliceomatic');

  % Context Menus
  % Slice Context Menu
  d.uic=uicontextmenu('callback', @slicecontextmenu);
  d.vistog = uimenu(d.uic,'label','Visible','callback','sliceomatic togglevisible');
  d.uicdelete = uimenu(d.uic,'label','Delete','callback','sliceomatic deleteslice');
  d.smfacet   = uimenu(d.uic,'label','Color Faceted','callback','sliceomatic setfaceted','separator','on');
  d.smflat    = uimenu(d.uic,'label','Color Flat','callback','sliceomatic setflat');
  d.sminterp  = uimenu(d.uic,'label','Color Interp','callback','sliceomatic setinterp');
  d.smtex     = uimenu(d.uic,'label','Color Texture','callback','sliceomatic settexture');
  d.smnone    = uimenu(d.uic,'label','Color None','callback','sliceomatic setnone');
  d.smtnone   = uimenu(d.uic,'label','Transparency None','callback','sliceomatic setalphanone','separator','on');
  d.smtp5     = uimenu(d.uic,'label','Transparency .5','callback','sliceomatic setalphapoint5');
  d.smtflat   = uimenu(d.uic,'label','Transparency Flat','callback','sliceomatic setalphaflat');
  d.smtinterp = uimenu(d.uic,'label','Transparency Interp','callback','sliceomatic setalphainterp');
  d.smttex    = uimenu(d.uic,'label','Transparency Texture','callback','sliceomatic setalphatexture');
  d.smcontour = uimenu(d.uic,'label','Add Contour','separator','on');
  d.smcont0   = uimenu(d.smcontour,'label','Auto','callback','sliceomatic slicecontour');
  d.smcont1   = uimenu(d.smcontour,'label','Select','callback','sliceomatic slicecontour_select','separator','on');
  d.smrcontour= uimenu(d.uic,'label','Remove Contour','callback','sliceomatic deleteslicecontour');
  d.smcflat   = uimenu(d.uic,'label','Contour Flat','callback','sliceomatic slicecontourflat','separator','on');
  d.smcinterp = uimenu(d.uic,'label','Contour Interp','callback','sliceomatic slicecontourinterp');
  d.smcblack  = uimenu(d.uic,'label','Contour Black','callback','sliceomatic slicecontourblack');
  d.smcwhite  = uimenu(d.uic,'label','Contour White','callback','sliceomatic slicecontourwhite');
  d.smccolor  = uimenu(d.uic,'label','Contour Color','callback','sliceomatic slicecontourcolor');
  d.smclinew  = uimenu(d.uic,'label','Contour Line Width');
  d.smcl1     = uimenu(d.smclinew,'label','1','callback','sliceomatic slicecontourlinewidth 1');
  d.smcl2     = uimenu(d.smclinew,'label','2','callback','sliceomatic slicecontourlinewidth 2');
  d.smcl3     = uimenu(d.smclinew,'label','3','callback','sliceomatic slicecontourlinewidth 3');
  d.smcl4     = uimenu(d.smclinew,'label','4','callback','sliceomatic slicecontourlinewidth 4');
  d.smcl5     = uimenu(d.smclinew,'label','5','callback','sliceomatic slicecontourlinewidth 5');
  d.smcl6     = uimenu(d.smclinew,'label','6','callback','sliceomatic slicecontourlinewidth 6');

  % Isosurface Context Menu
  d.uiciso=uicontextmenu('callback',@isocontextmenu);
  d.vistogiso = uimenu(d.uiciso,'label','Visible','callback','sliceomatic isotogglevisible');
  d.isodelete = uimenu(d.uiciso,'label','Delete','callback','sliceomatic isodelete');
  d.isoflatlight=uimenu(d.uiciso,'label','Lighting Flat','callback','sliceomatic isoflatlight','separator','on');
  d.isosmoothlight=uimenu(d.uiciso,'label','Lighting Smooth','callback','sliceomatic isosmoothlight');
  d.isocolor = uimenu(d.uiciso,'label','Change Color','callback','sliceomatic isocolor','separator','on');
  d.isoalpha=uimenu(d.uiciso,'label','Change Transparency');
  uimenu(d.isoalpha,'label','.2','callback','sliceomatic isoalpha .2');
  uimenu(d.isoalpha,'label','.5','callback','sliceomatic isoalpha .5');
  uimenu(d.isoalpha,'label','.8','callback','sliceomatic isoalpha .8');
  uimenu(d.isoalpha,'label','1','callback','sliceomatic isoalpha 1');
  d.isocap=uimenu(d.uiciso,'label','Add IsoCaps','callback','sliceomatic isocaps','separator','on');

  outd = d;

function controlmenu(fig, action)
% Handle doing things to the CONTROLS menu

  d=getappdata(gcf,'sliceomatic');

  if cameratoolbar('getvisible')
    set(d.camtoolbar,'checked','on');
  else
    set(d.camtoolbar,'checked','off');
  end

  set([d.dcalpha1 d.dcalpha8 d.dcalpha6 d.dcalpha5 d.dcalpha6 d.dcalpha2 d.dcalpha0...
       d.dclabels d.dcvis ],...
      'checked','off');

  switch get(d.pxx,'facealpha')
   case 1,  set(d.dcalpha1,'checked','on');
   case .8, set(d.dcalpha8,'checked','on');
   case .6, set(d.dcalpha6,'checked','on');
   case .5, set(d.dcalpha5,'checked','on');
   case .4, set(d.dcalpha4,'checked','on');
   case .2, set(d.dcalpha2,'checked','on');
   case 0,  set(d.dcalpha0,'checked','on');
  end

  if ~isempty(get(d.axx,'xticklabel'))
    set(d.dclabels,'checked','on');
  end

  if strcmp(get(d.axx,'visible'),'on')
    set(d.dcvis,'checked','on');
  end

  if 0
    xt = get(get(d.axx,'title'),'string');
    switch xt
     case 'X Slice Controller'
      set(d.dcslice,'checked','on');
    end

    xt = get(get(d.axiso,'title'),'string');
    switch xt
     case 'Iso Surface Controller'
      set(d.dciso,'checked','on');
    end
  end

function defaultmenu(fig, action)
% Handle toggling bits on the slice defaults menu

  d=getappdata(gcf,'sliceomatic');

  set([d.dfacet d.dflat d.dinterp d.dtex d.dtnone d.dtflat d.dtinterp ...
       d.dttex d.dcflat d.dcinterp d.dcblack d.dcwhite d.dcnone ...
       d.dlflat d.dlsmooth ...
       d.smcl1 d.smcl2 d.smcl3 d.smcl4 d.smcl5 d.smcl6 ], 'checked','off');
  switch d.defcolor
   case 'faceted'
    set(d.dfacet,'checked','on');
   case 'flat'
    set(d.dflat,'checked','on');
   case 'interp'
    set(d.dinterp,'checked','on');
   case 'texture'
    set(d.dtex,'checked','on');
   case 'none'
    set(d.dcnone,'checked','on');
  end
  switch d.defalpha
   case 'none'
    set(d.dtnone,'checked','on');
   case 'flat'
    set(d.dtflat,'checked','on');
   case 'interp'
    set(d.dtinterp,'checked','on');
   case 'texture'
    set(d.dttex,'checked','on');
  end
  switch d.deflight
   case 'flat'
    set(d.dlflat,'checked','on');
   case 'smooth'
    set(d.dlsmooth,'checked','on');
  end
  switch d.defcontourcolor
   case 'flat'
    set(d.dcflat,'checked','on');
   case 'interp'
    set(d.dcinterp,'checked','on');
   case 'black'
    set(d.dcblack,'checked','on');
   case 'white'
    set(d.dcwhite,'checked','on');
  end
  switch d.defcontourlinewidth
   case 1, set(d.dcl1,'checked','on');
   case 2, set(d.dcl2,'checked','on');
   case 3, set(d.dcl3,'checked','on');
   case 4, set(d.dcl4,'checked','on');
   case 5, set(d.dcl5,'checked','on');
   case 6, set(d.dcl6,'checked','on');
  end

function slicecontextmenu(fig,action)
% Context menu state for slices

  d=getappdata(gcf,'sliceomatic');

  [a s]=getarrowslice;
  set([d.smfacet d.smflat d.sminterp d.smtex d.smtnone d.smtp5 ...
       d.smtflat d.smtinterp d.smttex d.smnone],'checked','off');
  set(d.vistog,'checked',get(s,'visible'));

  if propcheck(s,'edgec',[0 0 0])
    set(d.smfacet,'checked','on');
  elseif propcheck(s,'facec','flat')
    set(d.smflat,'checked','on');
  end
  if propcheck(s,'facec','interp')
    set(d.sminterp,'checked','on');
  end
  if propcheck(s,'facec','texturemap')
    set(d.smtex,'checked','on');
  end
  if propcheck(s,'facec','none')
    set(d.smnone,'checked','on');
  end
  if propcheck(s,'facea',1)
    set(d.smtnone,'checked','on');
  end
  if propcheck(s,'facea',.5)
    set(d.smtp5,'checked','on');
  end
  if propcheck(s,'facea','flat')
    set(d.smtflat,'checked','on');
  end
  if propcheck(s,'facea','interp')
    set(d.smtinterp,'checked','on');
  end
  if propcheck(s,'facea','texturemap')
    set(d.smttex,'checked','on');
  end
  cm = [d.smcflat d.smcinterp d.smcblack d.smcwhite d.smccolor ...
       d.smcl1 d.smcl2 d.smcl3 d.smcl4 d.smcl5 d.smcl6 ];
  set(cm,'checked','off');
  if isempty(getappdata(s,'contour'))
    set(d.smcontour,'enable','on');
    set(d.smrcontour,'enable','off');
    set(cm,'enable','off');
  else
    set(d.smcontour,'enable','off')
    set(d.smrcontour,'enable','on')
    set(cm,'enable','on')
    c = getappdata(s,'contour');
    ec = get(c,'edgecolor');
    if isa(ec,'char')
      switch ec
       case 'flat'
        set(d.smcflat,'checked','on');
       case 'interp'
        set(d.smcinterp,'checked','on');
      end
    else
      if ec == [ 1 1 1 ]
        set(d.smcwhite,'checked','on');
      elseif ec == [ 0 0 0 ]
        set(d.smcblack,'checked','on');
      else
        set(d.smccolor,'checked','on');
      end
    end
    clw = get(c,'linewidth');
    switch clw
     case 1, set(d.smcl1,'checked','on');
     case 2, set(d.smcl2,'checked','on');
     case 3, set(d.smcl3,'checked','on');
     case 4, set(d.smcl4,'checked','on');
     case 5, set(d.smcl5,'checked','on');
     case 6, set(d.smcl6,'checked','on');
    end
  end

function isocontextmenu(fig,action)
% Context menu state for isosurfaces

  d=getappdata(gcf,'sliceomatic');

  [a s]=getarrowslice;
  if propcheck(s,'facelighting','flat')
    set(d.isoflatlight,'checked','on');
    set(d.isosmoothlight,'checked','off');
  else
    set(d.isoflatlight,'checked','off');
    set(d.isosmoothlight,'checked','on');
  end
  set(d.vistogiso,'checked',get(s,'visible'));
  if ~isempty(getappdata(s,'isosurfacecap'))
    set(d.isocap,'checked','on');
  else
    set(d.isocap,'checked','off');
  end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ Pxx, Pyy, Pxy, coh, pha, freq ] = func_coherence( f1, f2, nfft, Fs, filt, n_overlap )

%
% Computing power spectrum, cross spectrum, coherence and phase
%
% All input parameters are equivalent to csd or the other spectrum
% related function in MATLAB.
%
% This program use csd.m in Signal Processing Toolbox.
%
% Input
%	f1 and f2	input data
%	nfft		number of data for FFT
%	Fs		sampling frequency
%	filt		filter vector, hanning(nfft/2)
%	n_overlap	number of overlap for smoothing
%
% Output
%	Pxx		spectrum of f1
%	Pyy		spectrum of f2
%	Pxy		cross spectrum between f1 and f2
%	coh		coherence
%	pha		phase
%	freq		frequency vector
%
%======================================================================
% Terms:
%
%       Distributed under the terms of the terms of the BSD License
%
% Copyright:
%
%       Nobuhito Mori
%           Disaster Prevention Research Institue
%           Kyoto University
%           mori@oceanwave.jp
%
%========================================================================

[Pxx,freq] = csd( f1, f1, nfft, Fs, filt, n_overlap );
[Pyy,freq] = csd( f2, f2, nfft, Fs, filt, n_overlap );
[Pxy,freq] = csd( f1, f2, nfft, Fs, filt, n_overlap );

Kxy  = real( Pxy );
Qxy  = imag( Pxy );
coh  = Pxy.*conj(Pxy)./(Pxx.*Pyy);
pha  = atan2( Qxy, Kxy );


text:
%% PLANETOID CLASS (planetoid.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The Moon class is an iteration of the spheriod obstacle class aimed
% mostly providing a reference for satelite simulation.

% Author: James A. Douthwaite 09/02/2019

classdef planetoid < obstacle_spheroid
    properties
        inclination;    % Inclination angle
        orbit;          % Orbit
        orbitalSpeed;   % Oribital speed
        mass;           % Planetoid mass
        axialRate = 1;  % Rotational rate about its vertical axes(rad/s)
    end
    %% ///////////////////////// MAIN METHODS /////////////////////////////
    methods
        % Constructor
        function this = planetoid(varargin)
            % This function constructs the cuboid obstacle. The object must
            % be imported and represented with a global position and
            % velocity as all other objects are in OMAS.

            % Call the super class
            this = this@obstacle_spheroid(varargin);

            % IMPORT THE OBJECT'S GEOMETRY IF IT EXISTS
            [this.GEOMETRY] = OMAS_graphics.scale(this.GEOMETRY,this.radius);  % Scale

            % //////////////// Check for user overrides ///////////////////
            [this] = this.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % Setup - X = [x;x_dot]' 3D STATE VECTOR
        function [this] = setup(this,localXYZVelocity,localXYZrotations)
            % The state initialiser must be called 'initialise_localState'
            % and instead calls the 'initialise_3DVelocities' function in
            % this case.
            [this] = this.setup_3DVelocities(localXYZVelocity,localXYZrotations);
            % ADD A ROTATIONAL RATE
            this.localState(12) = this.axialRate; % The earth rotates a constant rate about its z interial axis
        end
        % Main
        function [this] = main(this,TIME,varargin)
            % SIMPLE UPDATE OF LOCAL STATE
            dt = TIME.dt;
            X = this.localState(1:6,1);
            U = this.localState(7:12,1);

            % SIMPLY INTEGRATE THE RATES
            [dXdt] = this.SingleIntegratorDynamics(X,U);
            eulerState = this.localState;
            eulerState(1:6,1) = this.localState(1:6,1) + dt*dXdt;

            % UPDATE THE GLOBAL PROPERTIES
            [this] = this.GlobalUpdate_3DVelocities(dt,eulerState);
        end
    end
end

repo: OpenMAS


text:
function net = cnn_imagenet_init(varargin)
% CNN_IMAGENET_INIT  Initialize a standard CNN for ImageNet

opts.scale = 1 ;
opts.initBias = 0.1 ;
opts.weightDecay = 1 ;
%opts.weightInitMethod = 'xavierimproved' ;
opts.weightInitMethod = 'gaussian' ;
opts.model = 'alexnet' ;
opts.batchNormalization = false ;
opts = vl_argparse(opts, varargin) ;

% Define layers
switch opts.model
  case 'alexnet'
    net.normalization.imageSize = [227, 227, 3] ;
    net = alexnet(net, opts) ;
  case 'vgg-f'
    net.normalization.imageSize = [224, 224, 3] ;
    net = vgg_f(net, opts) ;
  case 'vgg-m'
    net.normalization.imageSize = [224, 224, 3] ;
    net = vgg_m(net, opts) ;
  case 'vgg-s'
    net.normalization.imageSize = [224, 224, 3] ;
    net = vgg_s(net, opts) ;
  case 'vgg-vd-16'
    net.normalization.imageSize = [224, 224, 3] ;
    net = vgg_vd(net, opts) ;
  case 'vgg-vd-19'
    net.normalization.imageSize = [224, 224, 3] ;
    net = vgg_vd(net, opts) ;
  otherwise
    error('Unknown model ''%s''', opts.model) ;
end

% final touches
switch lower(opts.weightInitMethod)
  case {'xavier', 'xavierimproved'}
    net.layers{end}.weights{1} = net.layers{end}.weights{1} / 10 ;
end
net.layers{end+1} = struct('type', 'softmaxloss', 'name', 'loss') ;

net.normalization.border = 256 - net.normalization.imageSize(1:2) ;
net.normalization.interpolation = 'bicubic' ;
net.normalization.averageImage = [] ;
net.normalization.keepAspect = true ;

% --------------------------------------------------------------------
function net = add_block(net, opts, id, h, w, in, out, stride, pad, init_bias)
% --------------------------------------------------------------------
info = vl_simplenn_display(net) ;
fc = (h == info.dataSize(1,end) && w == info.dataSize(2,end)) ;
if fc
  name = 'fc' ;
else
  name = 'conv' ;
end
net.layers{end+1} = struct('type', 'conv', 'name', sprintf('%s%s', name, id), ...
                           'weights', {{init_weight(opts, h, w, in, out, 'single'), zeros(out, 1, 'single')}}, ...
                           'stride', stride, ...
                           'pad', pad, ...
                           'learningRate', [1 2], ...
                           'weightDecay', [opts.weightDecay 0]) ;
if opts.batchNormalization
  net.layers{end+1} = struct('type', 'bnorm', 'name', sprintf('bn%d',id), ...
                             'weights', {{ones(out, 1, 'single'), zeros(out, 1, 'single')}}, ...
                             'learningRate', [2 1], ...
                             'weightDecay', [0 0]) ;
end
net.layers{end+1} = struct('type', 'relu', 'name', sprintf('relu%s',id)) ;

% -------------------------------------------------------------------------
function weights = init_weight(opts, h, w, in, out, type)
% -------------------------------------------------------------------------
% See K. He, X. Zhang, S. Ren, and J. Sun. Delving deep into
% rectifiers: Surpassing human-level performance on imagenet
% classification. CoRR, (arXiv:1502.01852v1), 2015.

switch lower(opts.weightInitMethod)
  case 'gaussian'
    sc = 0.01/opts.scale ;
    weights = randn(h, w, in, out, type)*sc;
  case 'xavier'
    sc = sqrt(3/(h*w*in)) ;
    weights = (rand(h, w, in, out, type)*2 - 1)*sc ;
  case 'xavierimproved'
    sc = sqrt(2/(h*w*out)) ;
    weights = randn(h, w, in, out, type)*sc ;
  otherwise
    error('Unknown weight initialization method''%s''', opts.weightInitMethod) ;
end

% --------------------------------------------------------------------
function net = add_norm(net, opts, id)
% --------------------------------------------------------------------
if ~opts.batchNormalization
  net.layers{end+1} = struct('type', 'normalize', ...
                             'name', sprintf('norm%s', id), ...
                             'param', [5 1 0.0001/5 0.75]) ;
end

% --------------------------------------------------------------------
function net = add_dropout(net, opts, id)
% --------------------------------------------------------------------
if ~opts.batchNormalization
  net.layers{end+1} = struct('type', 'dropout', ...
                             'name', sprintf('dropout%s', id), ...
                             'rate', 0.5) ;
end

% --------------------------------------------------------------------
function net = alexnet(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;

net = add_block(net, opts, '1', 11, 11, 3, 96, 4, 0) ;
net = add_norm(net, opts, '1') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '2', 5, 5, 48, 256, 1, 2) ;
net = add_norm(net, opts, '2') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '3', 3, 3, 256, 384, 1, 1) ;
net = add_block(net, opts, '4', 3, 3, 192, 384, 1, 1) ;
net = add_block(net, opts, '5', 3, 3, 192, 256, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '6', 6, 6, 256, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

% --------------------------------------------------------------------
function net = vgg_s(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;
net = add_block(net, opts, '1', 7, 7, 3, 96, 2, 0) ;
net = add_norm(net, opts, '1') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 3, ...
                           'pad', [0 2 0 2]) ;

net = add_block(net, opts, '2', 5, 5, 96, 256, 1, 0) ;
net = add_norm(net, opts, '2') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ;

net = add_block(net, opts, '3', 3, 3, 256, 512, 1, 1) ;
net = add_block(net, opts, '4', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '5', 3, 3, 512, 512, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 3, ...
                           'pad', [0 1 0 1]) ;

net = add_block(net, opts, '6', 6, 6, 512, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

% --------------------------------------------------------------------
function net = vgg_m(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;
net = add_block(net, opts, '1', 7, 7, 3, 96, 2, 0) ;
net = add_norm(net, opts, '1') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '2', 5, 5, 96, 256, 2, 1) ;
net = add_norm(net, opts, '2') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ;

net = add_block(net, opts, '3', 3, 3, 256, 512, 1, 1) ;
net = add_block(net, opts, '4', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '5', 3, 3, 512, 512, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '6', 6, 6, 512, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

% --------------------------------------------------------------------
function net = vgg_f(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;
net = add_block(net, opts, '1', 11, 11, 3, 64, 4, 0) ;
net = add_norm(net, opts, '1') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ;

net = add_block(net, opts, '2', 5, 5, 64, 256, 1, 2) ;
net = add_norm(net, opts, '2') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '3', 3, 3, 256, 256, 1, 1) ;
net = add_block(net, opts, '4', 3, 3, 256, 256, 1, 1) ;
net = add_block(net, opts, '5', 3, 3, 256, 256, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '6', 6, 6, 256, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

% --------------------------------------------------------------------
function net = vgg_vd(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;
net = add_block(net, opts, '1_1', 3, 3, 3, 64, 1, 1) ;
net = add_block(net, opts, '1_2', 3, 3, 64, 64, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '2_1', 3, 3, 64, 128, 1, 1) ;
net = add_block(net, opts, '2_2', 3, 3, 128, 128, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '3_1', 3, 3, 128, 256, 1, 1) ;
net = add_block(net, opts, '3_2', 3, 3, 256, 256, 1, 1) ;
net = add_block(net, opts, '3_3', 3, 3, 256, 256, 1, 1) ;
if strcmp(opts.model, 'vgg-vd-19')
  net = add_block(net, opts, '3_4', 3, 3, 256, 256, 1, 1) ;
end
net.layers{end+1} = struct('type', 'pool', 'name', 'pool3', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '4_1', 3, 3, 256, 512, 1, 1) ;
net = add_block(net, opts, '4_2', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '4_3', 3, 3, 512, 512, 1, 1) ;
if strcmp(opts.model, 'vgg-vd-19')
  net = add_block(net, opts, '4_4', 3, 3, 512, 512, 1, 1) ;
end
net.layers{end+1} = struct('type', 'pool', 'name', 'pool4', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '5_1', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '5_2', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '5_3', 3, 3, 512, 512, 1, 1) ;
if strcmp(opts.model, 'vgg-vd-19')
  net = add_block(net, opts, '5_4', 3, 3, 512, 512, 1, 1) ;
end
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '6', 7, 7, 512, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

author: huangzehao

repo: caffe-vdsr


text:
%  H1 line
%
%  ::
%
%
%  Args:
%
%  Returns:
%     :
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function r = revaltrig(zz, zj, fj, wj, form)
%REVALTRIG   Evaluate rational function in trigonometric barycentric form.
%   R = REVALTRIG(ZZ, ZJ, FJ, WJ, FORM) returns R (vector of floats), the
%   values of the trigonometric barycentric rational function of form FORM
%   with support points ZJ, function values FJ, and barycentric weights WJ
%   evaluated at the points ZZ.
%
% See also AAATRIG, PRZTRIG, REVAL.

% Copyright 2018 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.
m = numel(zj);
zv = zz(:);                             % vectorize zz if necessary

% Define basis functions
if strcmp(form,'even')
    cst = @(x) cot(x);
elseif strcmp(form,'odd')
    cst = @(x) csc(x);
end

CC = cst(bsxfun(@minus, zv, zj.')/2); % Cauchy matrix
r = (CC*(wj.*fj))./(CC*wj);             % vector of values

if strcmp(form,'odd')
r(isinf(real(zv/1i)) & real(zv/1i)>0) = sum(wj.*fj.*exp(-1i*zj/2))./sum(wj.*exp(-1i*zj/2));
r(isinf(real(zv/1i)) & real(zv/1i)<0) = sum(wj.*fj.*exp(+1i*zj/2))./sum(wj.*exp(+1i*zj/2));
elseif strcmp(form,'even')
r(isinf(real(zv/1i))) = sum(wj.*fj)./sum(wj);
end

% Deal with NaN:
ii = find(isnan(r));
for jj = 1:length(ii)
    if ( isnan(zv(ii(jj))) || ~any(zv(ii(jj)) == zj) )
        % r(NaN) = NaN is fine.
        % The second case may happen if r(zv(ii)) = 0/0 at some point.
    else
        % Clean up values NaN = inf/inf at support points.
        % Find the corresponding node and set entry to correct value:
        r(ii(jj)) = fj(zv(ii(jj)) == zj);
    end
end

% Reshape to input format:
r = reshape(r, size(zz));

end % End of REVALTRIG().

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function c=comp_wfbt(f,wtNodes,rangeLoc,rangeOut,ext)
%COMP_WFBT Compute Wavelet Filterbank Tree
%   Usage:  c=comp_wfbt(f,wtNodes,rangeLoc,rangeOut,ext);
%
%   Input parameters:
%         f        : Input L*W array.
%         wtNodes  : Filterbank tree nodes (elementary filterbanks) in
%                    BF order. Length *nodeNo* cell array of structures.
%         rangeLoc : Idxs of each node terminal outputs. Length *nodeNo*
%                    cell array of vectors.
%         rangeOut : Output subband idxs of each node terminal outputs.
%         ext      : Type of the forward transform boundary handling.
%
%   Output parameters:
%         c        : Cell array of coefficients. Each element is one
%                    subband (matrix with W columns).
%

% Do non-expansve transform if ext=='per'
doPer = strcmp(ext,'per');
% Pre-allocated output
c = cell(sum(cellfun(@(rEl) numel(rEl),rangeOut)),1);

 ca = {f};
 % Go over all nodes in breadth-first order
 for jj=1:numel(wtNodes)
    % Load current filterbank
    wtNode = wtNodes{jj}.h(:);
    % Node filters to a cell array
    % hCell = cellfun(@(hEl) conj(flipud(hEl.h(:))),wtNode,'UniformOutput',0);
    hCell = cellfun(@(hEl) hEl.h(:),wtNode,'UniformOutput',0);
    % Node filters subs. factors
    a = wtNodes{jj}.a;
    % Node filters initial skips
    if(doPer)
       %offset = cellfun(@(hEl) 1-numel(hEl.h)-hEl.offset,wtNode);
       offset = cellfun(@(hEl) hEl.offset,wtNode);
    else
       offset = -(a-1);
    end

    % Run filterbank
    catmp=comp_filterbank_td(ca{1},hCell,a,offset,ext);
    % Pick what goes directy to the output...
    c(rangeOut{jj}) = catmp(rangeLoc{jj});
    % and save the rest.
    diffRange = 1:numel(hCell);
    diffRange(rangeLoc{jj}) = [];
    ca = [ca(2:end);catmp(diffRange)];
 end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [TW] = hpb2TW(hpb)
% Convert power from boiler horsepower to terawatts.
% Chad A. Greene 2012
TW = hpb*9809.5e-12;


text:
function [vals, pos] = minandmax(f)
%MINANDMAX   Global minimum and maximum of the SINGFUN F on [-1,1].
%   VALS = MINANDMAX(F) returns a 2-vector VALS = [MIN(F); MAX(F)] with the
%   global minimum and maximum of the SINGFUN F on [-1,1].
%
%   [VALS, POS] = MINANDMAX(F) returns also the 2-vector POS where the minimum
%   and maximum of F occur.
%
% See also MIN, MAX.

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

tol = chebfunpref().blowupPrefs.exponentTol;

if ( ~any(f.exponents) || all(abs(f.exponents) < tol) )

    % The function is actually smooth!
    [vals, pos] = minandmax(f.smoothPart);

else

    % Initialise:
    minF = [];
    maxF = [];
    minLoc = [];
    maxLoc = [];

    % Look for blow up at the left:
    if ( f.exponents(1) < -tol ) % Singularity at the left end.
        fVal = feval(f, -1);
        if ( fVal == inf )
            maxF = inf;
            maxLoc = -1;
        elseif ( fVal == -inf )
            minF = -inf;
            minLoc = -1;
        else
            % NaNs may occur and then we can not conclude anything.
            error('CHEBFUN:SINGFUN:minandmax:boundedness', ...
            ['Function has a singularity but isn''t infinite at the left ' ...
             'endpoint']);
        end
    end

    % Look for blow up at the right:
    if ( f.exponents(2) < -tol ) % Singularity at the right end.
        fVal = feval(f, 1);
        if ( fVal == inf )
            maxF = inf;
            maxLoc = 1;
        elseif ( fVal == -inf )
            minF = -inf;
            minLoc = 1;
        else
            % NaNs may occur and then we can not conclude anything.
            error('CHEBFUN:SINGFUN:minandmax:boundedness', ...
            ['Function has a singularity but isn''t infinite at the right ' ...
             'endpoint']);
        end
    end

    % If min or max is empty, then we need to do more work:
    if ( isempty(minF) || isempty(maxF) )
        % Find the roots of the derivative for local minima:
        fp = diff(f);
        r = roots(fp);
        % Append the end points and remove duplicates:
        r = unique([-1 ; r ; 1]);
        if ( isempty(maxF) )
            % Take the maximum of the local maxima:
            [maxF, maxIndex] = max(feval(f, r));
            maxLoc = r(maxIndex);
        end
        if ( isempty(minF) )
            % Take the minimum of the local minima:
            [minF, minIndex] = min(feval(f, r));
            minLoc = r(minIndex);
        end
    end
    vals = [ minF; maxF ];
    pos = [ minLoc; maxLoc ];

end

end


text:
function  dtiFindFGOI

%This functions takes [fgA, A, embbasis] (this is atlas), a set of nvecs in that
%atlas space that define an fgOI, and a new fibergroup fgS. The new fibergroup
%get embedded into the space, and those fibers that fall within the
%boundaries of the fgOI, are highlighted.

nvec=size(fgOI, 2);

%Compute embedding vectors for S;
E=dtiNewDataOntoEmbeddingVectors(A, S, embbasis, nvec);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function symb_pvec = polyprint(pvec)
%POLYPRINT Pretty print polynomial expression
%
% POLYPRINT is obsolete. Use SDISPLAY instead.

for pi = 1:size(pvec,1)
    for pj = 1:size(pvec,2)
        p = pvec(pi,pj);

        if isa(p,'double')
            symb_p = num2str(p);
        else
            LinearVariables = depends(p);
            x = recover(LinearVariables);
            exponent_p = full(exponents(p,x));
            names = cell(length(x),1);
            W = evalin('caller','whos');
            for i = 1:size(W,1)
                if strcmp(W(i).class,'sdpvar')% | strcmp(W(i).class,'lmi')
                    thevars = evalin('caller',W(i).name)    ;
                    if is(thevars,'scalar') & is(thevars,'linear') & length(getvariables(thevars))==1
                        index_in_p = find(ismember(LinearVariables,getvariables(thevars)));
                        if ~isempty(index_in_p)
                            names{index_in_p}=W(i).name;
                        end
                    end
                end
            end

            symb_p = '';
            if all(exponent_p(1,:)==0)
                symb_p = num2str(full(getbasematrix(p,0)));
                exponent_p = exponent_p(2:end,:);
            end

            for i = 1:size(exponent_p,1)
                coeff = full(getbasematrixwithoutcheck(p,i));
                switch coeff
                    case 1
                        coeff='+';
                    case -1
                        coeff = '-';
                    otherwise
                        if isreal(coeff)
                        if coeff >0
                            coeff = ['+' num2str2(coeff)];
                        else
                            coeff=[num2str2(coeff)];
                        end
                        else
                            coeff = ['+' '(' num2str2(coeff) ')' ];
                        end
                end
                symb_p = [symb_p coeff symbmonom(names,exponent_p(i,:))];
            end
            if symb_p(1)=='+'
                symb_p = symb_p(2:end);
            end
        end
        symb_pvec{pi,pj} = symb_p;
    end
end

function s = symbmonom(names,monom)
s = '';
for j = 1:length(monom)
    if abs( monom(j))>0
        s = [s names{j}];
        if monom(j)~=1
            s = [s '^' num2str(monom(j))];
        end
    end
end

function s = num2str2(x)
        s = num2str(full(x));
        if isequal(s,'1')
            s = '';
        end
        if isequal(s,'-1')
            s = '-';
        end


text:
% ALIGNK  calculate kernel alignment
%
% [X]=ALIGNK(K1,K2) calculates alignment between kernels K1 and K2.

function [x]=alignk(ik,ok)

x=sum(sum( ik .* ok)) / sqrt(sum(sum(ik.^2)) * sum(sum(ok.^2)) );


text:
function [lm lm_avg] = rcu_ach(n, delta, epsil, plow, pup);
% Compute RCU achievability bound.
% lm returns the bound valid for maximal probability of error
% lm_avg is for the average
%
% This is not a vectorized version, because the function is really slow.
% You can speed it up significantly if you provide a bracket plow, pup for the value of logm.
% For example, a smart idea is to set
% plow = gallager_ach(n, delta, epsil) and pup = converse(n,delta,epsil);

if nargin < 5
	plow = 0;
	pup = n;
end

eps_test = precise_rand(n, delta, pup);
if(eps_test < epsil)
	disp(sprintf([	'-- achiev_prand(n = %d, delta = %g, epsil = %g): eps_test = %g\n'...
			'        This is a bug? precise_rand() contradicts converse ?!?!'], ...
			n, delta, epsil, eps_test));
	error('achiev_prand');
	lm_avg = 0; lm = 0;
	return;
end

eps_test = precise_rand(n, delta, plow);
if(eps_test > epsil)
	disp(sprintf([	'-- achiev_prand(n = %d, delta = %g, epsil = %g): eps_test = %g\n'...
			'        Can not find lower bound for precise_rand() !!!!'], ...
			n, delta, epsil, eps_test));
	%error('achiev_prand');
	lm_avg = 0; lm = 0;
	return;
end

% Take into account that we are computing AVERAGE prob. of error,
% and so we need a random linear code trick to go to maximum => log M must be integer
while floor(plow) < floor(pup);
	ptest = (plow + pup)/2;
	eps_test = precise_rand(n, delta, ptest);
	if(eps_test > epsil)
		pup = ptest;
	else
		plow = ptest;
		% This is not needed as we stop on floor-test
		%if  (epsil-eps_test)< 1e-2*epsil
		%	break;
		%end
	end
end
lm_avg = plow;
lm = floor(plow);

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function y = logSt(X, mu, sigma, v)
% Compute log pdf of a Student's t distribution.
% Input:
%   X: d x n data matrix
%   mu: mean
%   sigma: variance
%   v: degree of freedom
% Output:
%   y: probability density in logrithm scale y=log p(x)
% Written by mo Chen (sth4nth@gmail.com).
[d,k] = size(mu);

if size(sigma,1)==d && size(sigma,2)==d && k==1
    [R,p]= chol(sigma);
    if p ~= 0
        error('ERROR: sigma is not SPD.');
    end
    X = bsxfun(@minus,X,mu);
    Q = R'\X;
    q = dot(Q,Q,1);  % quadratic term (M distance)
    o = -log(1+q/v)*((v+d)/2);
    c = gammaln((v+d)/2)-gammaln(v/2)-(d*log(v*pi)+2*sum(log(diag(R))))/2;
    y = c+o;
elseif size(sigma,1)==d && size(sigma,2)==k
    lambda = 1./sigma;
    ml = mu.*lambda;
    q = bsxfun(@plus,X'.^2*lambda-2*X'*ml,dot(mu,ml,1)); % M distance
    o = bsxfun(@times,log(1+bsxfun(@times,q,1./v)),-(v+d)/2);
    c = gammaln((v+d)/2)-gammaln(v/2)-(d*log(pi*v)+sum(log(sigma),1))/2;
    y = bsxfun(@plus,o,c);
elseif size(sigma,1)==1 && size(sigma,2)==k
    X2 = repmat(dot(X,X,1)',1,k);
    D = bsxfun(@plus,X2-2*X'*mu,dot(mu,mu,1));
    q = bsxfun(@times,D,1./sigma);  % M distance
    o = bsxfun(@times,log(1+bsxfun(@times,q,1./v)),-(v+d)/2);
    c = gammaln((v+d)/2)-gammaln(v/2)-d*log(pi*v.*sigma)/2;
    y = bsxfun(@plus,o,c);
else
    error('Parameters are mismatched.');
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function stage_time_struct = timesStruct(N, nt)
stage_time_struct = ocl.types.Structure();
stage_time_struct.addRepeated({'states', 'integrator', 'controls'}, ...
  {ocl.types.Matrix([1,1]), ocl.types.Matrix([nt,1]), ocl.types.Matrix([1,1])}, N);
stage_time_struct.add('states', ocl.types.Matrix([1,1]));

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
% addpath(genpath('fm'));

m = 10;
S = zeros(m,m,m); S(2:4,2:4,2:4)=1;

options.null = 0;

if 0
options.constraint_map = S;
options.constraint_map(S==0) = -Inf;
options.constraint_map(S~=0) = +Inf;
end

% gaussian weight (path will avoid center of the cube)
x = -1:2/(m-1):1;
[X,Y,Z] = meshgrid(x,x,x);
sigma = 0.4;
W = 1./(1 + exp( -(X.^2+Y.^2+Z.^2)/sigma^2 ) );

% options.nb_iter_max = Inf;

W = rescale(S,1e-5,1);
W = rand(m,m,m);
W = rescale(W,.5,1);

tic
[D,RS] = perform_fast_marching(W, [2; 2; 2], options);
fprintf(1,'m = %d, time = %f\n', m, toc);

return;

for m=[10 20 50 100]

    S = zeros(m,m,m); S(2:4,2:4,2:4)=1;

    options.constraint_map = S;
    options.constraint_map(S==0) = -Inf;
    options.constraint_map(S~=0) = +Inf;

    W = rescale(S,1e-5,1);

    tic
    [D,RS] = perform_fast_marching(W, [2; 2; 2], options);
    fprintf(1,'m = %d, time = %f\n', m, toc);
end

author: alecjacobson

repo: gptoolbox


text:
function eps2pdf(source, dest, crop, append, gray, quality, gs_options)
%EPS2PDF  Convert an eps file to pdf format using ghostscript
%
% Examples:
%   eps2pdf source dest
%   eps2pdf(source, dest, crop)
%   eps2pdf(source, dest, crop, append)
%   eps2pdf(source, dest, crop, append, gray)
%   eps2pdf(source, dest, crop, append, gray, quality)
%   eps2pdf(source, dest, crop, append, gray, quality, gs_options)
%
% This function converts an eps file to pdf format. The output can be
% optionally cropped and also converted to grayscale. If the output pdf
% file already exists then the eps file can optionally be appended as a new
% page on the end of the eps file. The level of bitmap compression can also
% optionally be set.
%
% This function requires that you have ghostscript installed on your
% system. Ghostscript can be downloaded from: http://www.ghostscript.com
%
% Inputs:
%   source  - filename of the source eps file to convert. The filename is
%             assumed to already have the extension ".eps".
%   dest    - filename of the destination pdf file. The filename is assumed
%             to already have the extension ".pdf".
%   crop    - boolean indicating whether to crop the borders off the pdf.
%             Default: true.
%   append  - boolean indicating whether the eps should be appended to the
%             end of the pdf as a new page (if the pdf exists already).
%             Default: false.
%   gray    - boolean indicating whether the output pdf should be grayscale
%             or not. Default: false.
%   quality - scalar indicating the level of image bitmap quality to
%             output. A larger value gives a higher quality. quality > 100
%             gives lossless output. Default: ghostscript prepress default.
%   gs_options - optional ghostscript options (e.g.: '-dNoOutputFonts'). If
%                multiple options are needed, enclose in call array: {'-a','-b'}

% Copyright (C) Oliver Woodford 2009-2014, Yair Altman 2015-

% Suggestion of appending pdf files provided by Matt C at:
% http://www.mathworks.com/matlabcentral/fileexchange/23629

% Thank you to Fabio Viola for pointing out compression artifacts, leading
% to the quality setting.
% Thank you to Scott for pointing out the subsampling of very small images,
% which was fixed for lossless compression settings.

% 9/12/2011 Pass font path to ghostscript.
% 26/02/15: If temp dir is not writable, use the dest folder for temp
%           destination files (Javier Paredes)
% 28/02/15: Enable users to specify optional ghostscript options (issue #36)
% 01/03/15: Upon GS error, retry without the -sFONTPATH= option (this might solve
%           some /findfont errors according to James Rankin, FEX Comment 23/01/15)
% 23/06/15: Added extra debug info in case of ghostscript error; code indentation
% 04/10/15: Suggest a workaround for issue #41 (missing font path; thanks Mariia Fedotenkova)
% 22/02/16: Bug fix from latest release of this file (workaround for issue #41)
% 20/03/17: Added informational message in case of GS croak (issue #186)

    % Intialise the options string for ghostscript
    options = ['-q -dNOPAUSE -dBATCH -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress -sOutputFile="' dest '"'];
    % Set crop option
    if nargin < 3 || crop
        options = [options ' -dEPSCrop'];
    end
    % Set the font path
    fp = font_path();
    if ~isempty(fp)
        options = [options ' -sFONTPATH="' fp '"'];
    end
    % Set the grayscale option
    if nargin > 4 && gray
        options = [options ' -sColorConversionStrategy=Gray -dProcessColorModel=/DeviceGray'];
    end
    % Set the bitmap quality
    if nargin > 5 && ~isempty(quality)
        options = [options ' -dAutoFilterColorImages=false -dAutoFilterGrayImages=false'];
        if quality > 100
            options = [options ' -dColorImageFilter=/FlateEncode -dGrayImageFilter=/FlateEncode -c ".setpdfwrite << /ColorImageDownsampleThreshold 10 /GrayImageDownsampleThreshold 10 >> setdistillerparams"'];
        else
            options = [options ' -dColorImageFilter=/DCTEncode -dGrayImageFilter=/DCTEncode'];
            v = 1 + (quality < 80);
            quality = 1 - quality / 100;
            s = sprintf('<< /QFactor %.2f /Blend 1 /HSample [%d 1 1 %d] /VSample [%d 1 1 %d] >>', quality, v, v, v, v);
            options = sprintf('%s -c ".setpdfwrite << /ColorImageDict %s /GrayImageDict %s >> setdistillerparams"', options, s, s);
        end
    end
    % Enable users to specify optional ghostscript options (issue #36)
    if nargin > 6 && ~isempty(gs_options)
        if iscell(gs_options)
            gs_options = sprintf(' %s',gs_options{:});
        elseif ~ischar(gs_options)
            error('gs_options input argument must be a string or cell-array of strings');
        else
            gs_options = [' ' gs_options];
        end
        options = [options gs_options];
    end
    % Check if the output file exists
    if nargin > 3 && append && exist(dest, 'file') == 2
        % File exists - append current figure to the end
        tmp_nam = tempname;
        try
            % Ensure that the temp dir is writable (Javier Paredes 26/2/15)
            fid = fopen(tmp_nam,'w');
            fwrite(fid,1);
            fclose(fid);
            delete(tmp_nam);
        catch
            % Temp dir is not writable, so use the dest folder
            [dummy,fname,fext] = fileparts(tmp_nam); %#ok<ASGLU>
            fpath = fileparts(dest);
            tmp_nam = fullfile(fpath,[fname fext]);
        end
        % Copy the file
        copyfile(dest, tmp_nam);
        % Add the output file names
        options = [options ' -f "' tmp_nam '" "' source '"'];
        try
            % Convert to pdf using ghostscript
            [status, message] = ghostscript(options);
        catch me
            % Delete the intermediate file
            delete(tmp_nam);
            rethrow(me);
        end
        % Delete the intermediate file
        delete(tmp_nam);
    else
        % File doesn't exist or should be over-written
        % Add the output file names
        options = [options ' -f "' source '"'];
        % Convert to pdf using ghostscript
        [status, message] = ghostscript(options);
    end
    % Check for error
    if status
        % Retry without the -sFONTPATH= option (this might solve some GS
        % /findfont errors according to James Rankin, FEX Comment 23/01/15)
        orig_options = options;
        if ~isempty(fp)
            options = regexprep(options, ' -sFONTPATH=[^ ]+ ',' ');
            status = ghostscript(options);
            if ~status, return; end  % hurray! (no error)
        end
        % Report error
        if isempty(message)
            error('Unable to generate pdf. Check destination directory is writable.');
        elseif ~isempty(strfind(message,'/typecheck in /findfont'))
            % Suggest a workaround for issue #41 (missing font path)
            font_name = strtrim(regexprep(message,'.*Operand stack:\s*(.*)\s*Execution.*','$1'));
            fprintf(2, 'Ghostscript error: could not find the following font(s): %s\n', font_name);
            fpath = fileparts(mfilename('fullpath'));
            gs_fonts_file = fullfile(fpath, '.ignore', 'gs_font_path.txt');
            fprintf(2, '  try to add the font''s folder to your %s file\n\n', gs_fonts_file);
            error('export_fig error');
        else
            fprintf(2, '\nGhostscript error: perhaps %s is open by another application\n', dest);
            if ~isempty(gs_options)
                fprintf(2, '  or maybe the%s option(s) are not accepted by your GS version\n', gs_options);
            end
            fprintf(2, '  or maybe you have another gs executable in your system''s path\n');
            fprintf(2, 'Ghostscript options: %s\n\n', orig_options);
            error(message);
        end
    end
end

% Function to return (and create, where necessary) the font path
function fp = font_path()
    fp = user_string('gs_font_path');
    if ~isempty(fp)
        return
    end
    % Create the path
    % Start with the default path
    fp = getenv('GS_FONTPATH');
    % Add on the typical directories for a given OS
    if ispc
        if ~isempty(fp)
            fp = [fp ';'];
        end
        fp = [fp getenv('WINDIR') filesep 'Fonts'];
    else
        if ~isempty(fp)
            fp = [fp ':'];
        end
        fp = [fp '/usr/share/fonts:/usr/local/share/fonts:/usr/share/fonts/X11:/usr/local/share/fonts/X11:/usr/share/fonts/truetype:/usr/local/share/fonts/truetype'];
    end
    user_string('gs_font_path', fp);
end

author: optimaltransport


text:
function [flag] = point_visible(V,F,o,heuristic_ratio)
  % POINT_VISIBLE  test whether vertices of mesh are visible to a set of
  % points
  %
  % [flag] = point_visible(V,F,o)
  %
  % Input:
  %    V  #V by 3 list of vertex positions
  %    F  #F by 3 list of triangle indices
  %    o  row vector of position to test
  % Output:
  %    flag  #V by 1 list of bools (true) visible, (false) obstructed
  %

  dim = size(V,2);
  % Try to use accelerated bone_visible code
  if dim == 3 && 3==exist('bone_visible_embree','file')
    flag = bone_visible_embree(V,F,o,o);
  % otherwise try to use mex bone_visible code
  elseif dim ==3 && 3==exist('bone_visible','file')
    warning('Using non accelerated visibility test');
    flag = bone_visible(V,F,o,o);
  % otherwise use super slow matlab code
  else
    warning('Using non accelerated, pure matlab visibility test');

    assert(size(o,2) == dim);
    assert(size(F,2) == 3);
    % number of mesh vertices
    nv =  size(V,1);

    if dim == 2
      % get polygon edges of outline of 2D mesh
      O = outline(F);
    end

    if(exist('heuristic_ratio','var'))
      indices = randperm(nv);
      if(heuristic_ratio < 1)
        indices = indices(1:round(heuristic_ratio*nv));
      else
        indices = indices(1:round(heuristic_ratio));
      end
      flag = 0.5+eps + zeros(nv,1);
    else
      indices = 1:nv;
      flag = zeros(nv,1);
    end

    % loop over mesh vertices
    for jj_i = 1:numel(indices)
      progressbar(jj_i,numel(indices));
      jj = indices(jj_i);
      % query point, can o see q ?
      vjj = V(jj,:);

      dir = vjj-o;

      if dim == 3
        % extract triangles not containing jj
        Fmjj = F( all(F~=jj,2),:);
        %Fmjj = F;
        % shoot ray from o to q and see if any triangles not containing query get
        % hit
        %[hitc,tc] = mexray_mesh_intersect(o,dir,V,Fmjj);
        [hit,t] = ray_mesh_intersect(o,dir,V,Fmjj);
      elseif dim == 2
        Omjj = O( all(O~=jj,2),:);
        [hit,t] = ray_polygon_intersect(o,dir,V,Omjj);
      else
        error('Bad dimension');
      end

      %vjj-o
      %meshplot(V,F,'FC',repmat(hit,1,3))
      dir_mag =  sqrt(sum(dir.^2,2));
      flag(jj) = ~any(t(hit) < dir_mag);
    end

    if(exist('heuristic_ratio','var'))
      L = adjacency_matrix(F);
      L = L - diag(sum(L));
      i = 0;
      while(any(flag == (0.5+eps)))
        flag = flag+0.08*L*flag;
        flag(flag ~= (0.5+eps)) = round(flag(flag ~= (0.5+eps)));
        i = i +1;
      end
      i
      %L = cotmatrix(V,F);
      %I = speye(nv,nv);
      %L(indices,:) = I(indices,:);
      %flag = (L\flag)>0.5;
    end

  end
end


text:
function [fc,cc,kc,alpha_c,Rc,Tc,omc,nx,ny,x_dist,xd] = willson_convert(Ncx,Nfx,dx,dy,dpx,dpy,Cx,Cy,sx,f,kappa1,Tx,Ty,Tz,Rx,Ry,Rz,p1,p2);

%Conversion from Reg Willson's calibration format to my format

% Conversion:

% Focal length:
fc = [sx/dpx ; 1/dpy]*f;

% Principal point;
cc = [Cx;Cy];

% Skew:
alpha_c = 0;

% Extrinsic parameters:
Rx = rodrigues([Rx;0;0]);
Ry = rodrigues([0;Ry;0]);
Rz = rodrigues([0;0;Rz]);

Rc = Rz * Ry * Rx;

omc = rodrigues(Rc);

Tc = [Tx;Ty;Tz];

% More tricky: Take care of the distorsion:

Nfy = round(Nfx * 3/4);

nx = Nfx;
ny = Nfy;

% Select a set of DISTORTED coordinates uniformely distributed across the image:

[xp_dist,yp_dist] = meshgrid(0:Nfx-1,0:Nfy);

xp_dist = xp_dist(:)';
yp_dist = yp_dist(:)';

% Apply UNDISTORTION according to Willson:

xp_sensor_dist = dpx*(xp_dist - Cx)/sx;
yp_sensor_dist = dpy*(yp_dist - Cy);

dist_fact = 1 + kappa1*(xp_sensor_dist.^2 + yp_sensor_dist.^2);

xp_sensor = xp_sensor_dist .* dist_fact;
yp_sensor = yp_sensor_dist .* dist_fact;

xp = xp_sensor * sx / dpx  + Cx;
yp = yp_sensor / dpy  + Cy;

ind= find((xp > 0) & (xp < Nfx-1) & (yp > 0) & (yp < Nfy-1));

xp = xp(ind);
yp = yp(ind);
xp_dist = xp_dist(ind);
yp_dist = yp_dist(ind);

% Now, find my own set of parameters:

x_dist = [(xp_dist - cc(1))/fc(1);(yp_dist - cc(2))/fc(2)];
x_dist(1,:) = x_dist(1,:) - alpha_c * x_dist(2,:);

x = [(xp - cc(1))/fc(1);(yp - cc(2))/fc(2)];
x(1,:) = x(1,:) - alpha_c * x(2,:);

k = [0;0;0;0;0];

for kk = 1:5,

	[xd,dxddk] = apply_distortion(x,k);

	err = x_dist - xd;

	%norm(err)

   k_step = inv(dxddk'*dxddk)*(dxddk')*err(:);

   k = k + k_step; %inv(dxddk'*dxddk)*(dxddk')*err(:);

   %norm(k_step)/norm(k)

   if norm(k_step)/norm(k) < 10e-10,
      break;
   end;

end;

kc = k;

repo: ekfmonoslam