text:
randn('seed', 1e6);
rand('seed', 1e6);
if ~exist('mappingKern', 'var'), mappingKern = 'rbfardjit'; end
if ~exist('indPoints', 'var'), indPoints = 20; end
if ~exist('latentDim', 'var'), latentDim = 10; end
% Define a temporal model
if ~exist('dynamicsConstrainType', 'var'), dynamicsConstrainType = {'time'}; end
if ~exist('initVardistIters', 'var'), initVardistIters = 500; end
if ~exist('itNo', 'var'), itNo = [50 50]; end
if ~exist('dynamicKern', 'var'), dynamicKern = {'matern32','white','bias'}; end

%%
% load data
[Y, lbls] = lvmLoadData('brendan');
Y = Y(801:880,:);
t = linspace(1,2*pi, size(Y,1))';

h = 20;
w = 28;
scrsz = get(0,'ScreenSize');
figure('Position',[0.2*scrsz(3) 0.2*scrsz(4) 20*5 28*5])
for i=1:size(Y,1)
    imagesc(reshape(Y(i,:), h, w)'); colormap('gray');
    pause(0.05);
end

%% Split between training and test blocks
mask = [];
lastTrPts = 5;
r=1; % start with tr. set
while length(mask)<size(Y,1)-lastTrPts %The last lastTrPts will be from YTr necessarily
    blockSize = randperm(8);
    blockSize = blockSize(1);
    pts = min(blockSize, size(Y,1)-lastTrPts - length(mask));
    if r
        mask = [mask ones(1,pts)];
    else
        mask = [mask zeros(1,pts)];
    end
    r = ~r; % alternate between tr. and test set
end
mask = [mask ones(1,lastTrPts)];
indTr = find(mask);
indTs = find(~mask);
if sum(sort([indTr indTs]) - (1:size(Y,1)))
    error('Something went wrong in the dataset splitting...');
end
Nstar = length(indTs);

Ytr = Y(indTr,:); t_tr = t(indTr, :);
Yts = Y(indTs,:); t_ts = t(indTs, :);

% Additionally define some pixels to be observed for the test set
cutPoint=round(h/2)+1;
mask = zeros(h,1);
mask(cutPoint) = 1;
mask = [ones(1,cutPoint) zeros(1,h-cutPoint)];
mask=repmat(mask, 1,w);
indexMissing = find(mask);

YtsOriginal = Yts;
Yts(:, indexMissing) = NaN;
% Play the test set fully observed (left) and partially observed (right)
scrsz = get(0,'ScreenSize');
figure('Position',[0.2*scrsz(3) 0.2*scrsz(4) 20*15 28*5])
for i=1:size(Yts,1)
    subplot(1,2,1)
    imagesc(reshape(YtsOriginal(i,:), h, w)'); colormap('gray');
    subplot(1,2,2);
    imagesc(reshape(Yts(i,:), h, w)'); colormap('gray');
    pause(0.05);
end
%% Train models

% Options for the model
vargplvm_init; % Returns a configuration structure 'globalOpt'
options = vargplvmOptions('dtcvar');
options.kern = mappingKern;
options.numActive = indPoints;
options.optimiser = 'scg2';
options.latentDim = latentDim;
options.initSNR = 100;
if ~isempty('dynamicsConstrainType')
    % Temporal model (VGPDS)
    optionsDyn.type = 'vargpTime';
    optionsDyn.t=t;
    if ~isstruct(dynamicKern)
        optionsDyn.kern = kernCreate(t, dynamicKern);
    else
        optionsDyn.kern = dynamicKern;
    end
    % Create and optimise the model
    [~, ~, ~, ~, modelInitVardist] = vargplvmEmbed2(Ytr, latentDim, options, initVardistIters, 0, true, optionsDyn);
     model = vargplvmOptimiseModel(modelInitVardist, true, true, {0,itNo}, true);
else
    % Non-dynamical model (Bayesian GP-LVM)
    % Create and optimise the model
    [~, ~, ~, ~, modelInitVardist] = vargplvmEmbed2(Ytr, latentDim, options, initVardistIters, 0, true);
    model = vargplvmOptimiseModel(modelInitVardist, true, true, {0,itNo}, true);
end

%%

fprintf(1, '# Predicting only with the test time points...\n')
[Testmeans2 Testcovars2] = vargplvmPredictPoint(model.dynamics, t_ts);
Varmu2 = vargplvmPosteriorMeanVar(model, Testmeans2, Testcovars2);
% Mean absolute error per pixel
errorOnlyTimes = mean(abs(Varmu2(:) - YtsOriginal(:)));

%%
% Play the test set fully observed (left) and predicted (right)
scrsz = get(0,'ScreenSize');
figure('Position',[0.2*scrsz(3) 0.2*scrsz(4) 20*15 28*5])
t_ind = 1:size(Y,1);
t_ind = t_ind(indTs);
for i=1:size(Yts,1)
    subplot(1,2,1)
    imagesc(reshape(YtsOriginal(i,:), h, w)'); colormap('gray'); title(num2str(t_ind(i)))
    subplot(1,2,2);
    imagesc(reshape(Varmu2(i,:), h, w)'); colormap('gray');
    pause;
end


text:
function noise = ncnmNoiseParamInit(noise, y)

% NCNMNOISEPARAMINIT null category noise model's parameter initialisation.
% The null category noise model enables semi-supervised learning
% with Gaussian processes. The approach is described in a 2004 NIPS
% paper by Lawrence and Jordan.
%
% FORMAT
% DESC initialises the parameters of the null category noise model.
% ARG noise : the structure to initialise.
% ARG y : a set of target values.
% RETURN noise : the initialised noise structure.
%
% FORMAT
% DESC initialises the parameters of the null category noise model.
% ARG noise : the structure to initialise.
% RETURN noise : the initialised noise structure.
%
% SEEALSO : noiseParamInit, noiseCreate
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% NOISE

% The likelihood is not log concave.
noise.logconcave = 0;
noise.gammaSplit = 0;

if nargin > 1
  nClass1 = sum(y==1, 1);
  nClass2 = sum(y==-1, 1);
  totClass = nClass1 + nClass2;
  p1 = nClass1./totClass;
  noise.numProcess = size(y, 2);
  noise.gamman = sum(isnan(y))/length(y);
  noise.gammap = noise.gamman;
  noise.bias = invCumGaussian(p1);
else
  noise.bias = zeros(1, noise.numProcess);
  noise.gamman = 0.5;
  noise.gammap = 0.5;
end
if noise.gammaSplit
  noise.nParams = noise.numProcess+2;
else
  noise.nParams = noise.numProcess+1;
end

% Constrain noise.prior to be between 0 and 1.
if noise.gammaSplit
  noise.transforms.index = [noise.numProcess+1 noise.numProcess+2];
else
  noise.transforms.index = [noise.numProcess+1];
end
noise.transforms.type = optimiDefaultConstraint('zeroone');

% This isn't optimised, it sets the gradient of the erf.
noise.sigma2 = eps;

% Can handle missing values?
noise.missing = 1;
noise.width = 1;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% Import Script for EBSD Data
%
% This script was automatically created by the import wizard. You should
% run the whoole script or parts of it in order to import your data. There
% is no problem in making any changes to this script.

%% Specify Crystal and Specimen Symmetries

% crystal symmetry
CS = {crystal symmetry};

% plotting convention
setMTEXpref('xAxisDirection',{xAxisDirection});
setMTEXpref('zAxisDirection',{zAxisDirection});

%% Specify File Names

% path to files
pname = {path to files};

% which files to be imported
fname = {file names};

%% Z-Values

Z = {Z-values};

%% Import the Data

% create an EBSD variable containing the data
ebsd = EBSD.load(fname,CS,'interface',{interface},{Z},{options});

%% Correct Data

rot = rotation.byEuler({phi1},{Phi},{phi2});
ebsd = rotate(ebsd,rot,{rotationOption});

author: mtex-toolbox

repo: mtex

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


text:
%% ppval_extrapVal
% Below is a demonstration of the features of the |ppval_extrapVal| function

%%
clear; close all; clc;

%% Syntax
% |v=ppval_extrapVal(pp,xx,interpLim,extrapVal);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
        % Incarcarea tabelului
x=1:10;
y=1:2:19;
tab=[x',y']
        % Cautarea in tabel
a=table1(tab,2.5)

author: Sable

repo: mcbench-benchmarks


text:
function orig_path = mrtrix_set_ld_path
% Set the library path for mrtrix libraries

% Save origional environment path
orig_path = getenv('LD_LIBRARY_PATH');
% Get mrtrix path
[status, pathstr] = system('which csdeconv');
if status~=0
    error('Please install mrtrix')
end
% This will be the path to the mrtrix libraries
pathstr = fullfile(pathstr(1:end-13),'lib');
% set the environment path
setenv('LD_LIBRARY_PATH',pathstr)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function anglepairs = samplerandfeat(num_feats)
%SAMPLERANDFEAT Summary of this function goes here
%   Function: generate the locations of pixel pairs randomly
%   Detailed explanation goes here
%   Input:
%        num_feats: number of features
%        max_radius: the maximum radius of local region
%   Output:
%         anglepairs: the angles of pixel pairs

thetas_a = 2*pi*[0:1/(num_feats-1):1];
thetas_b = 2*pi*[0:1/(num_feats-1):1];

anglepairs = [thetas_a(randperm(length(thetas_a)))' thetas_b(randperm(length(thetas_b)))'];

end

author: jwyang

repo: face-alignment


text:
function Offspring = Operator(Problem,Particle,Pbest,Gbest)
% <operator> <real>
% Particle swarm optimization in NMPSO

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Parameter setting
    ParticleDec = Particle.decs;
    PbestDec    = Pbest.decs;
    GbestDec    = Gbest.decs;
    [N,D]       = size(ParticleDec);
    ParticleVel = Particle.adds(zeros(N,D));

    %% Particle swarm optimization
    W  = repmat(unifrnd(0.1,0.5,N,1),1,D);
    r1 = repmat(rand(N,1),1,D);
    r2 = repmat(rand(N,1),1,D);
    r3 = repmat(rand(N,1),1,D);
    C1 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    C2 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    C3 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    OffVel    = W.*ParticleVel + C1.*r1.*(PbestDec-ParticleDec) + C2.*r2.*(GbestDec-ParticleDec) + C3.*r3.*(GbestDec-PbestDec);
    OffDec    = ParticleDec + OffVel;
    Offspring = Problem.Evaluation(OffDec,OffVel);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Get vertex coordinates for the patch type 'B'
% [C, N, S, W, E] = HealpixGetPatchVertexCoordsB(n, i, j, INFO)
%
% Parameters
% n : grid resolution
% i : ring index
% j : intra-ring index
% INFO : intermediate information (output of HealpixSelectPatchClass())
% C : intra-patch coordinates
% N : coordinates for north vertex
% S : coordinates for south vertex
% W : coordinates for west vertex
% E : coordinates for east vertex

function [C, N, S, W, E] = HealpixGetPatchVertexCoordsB(n, i, j, INFO)

% gradient of the border
grad = INFO.polar_part;

decimal_i = INFO.decimal_i_n;
decimal_j = INFO.decimal_polar_intra_part;

int_i = INFO.int_i_n;

if int_i == n
    int_j = fix(j);
    if int_j == 0
        int_j = 4 * INFO.int_i_n;
    end

    north_int_j = mod(int_j - grad - 1, 4 * (n - 1)) + 1;
    east_int_j = mod(int_j + 1 - 1, 4 * n) + 1;

    N = [int_i - 1, north_int_j];
    S = [int_i + 1, int_j];
    W = [int_i    , int_j];
    E = [int_i    , east_int_j];

    offset_j = 0.5 * decimal_i;
    C = [decimal_i, decimal_j + offset_j - 0.5];
else
    int_j = fix(j - grad * decimal_i);
    if int_j == 0
        int_j = 4 * n;
    end

    west_int_j = mod(int_j + grad - 1, 4 * n) + 1;
    east_int_j = mod(west_int_j + 1 - 1, 4 * n) + 1;

    N = [int_i    , int_j];
    S = [int_i + 2, west_int_j];
    W = [int_i + 1, west_int_j];
    E = [int_i + 1, east_int_j];

    offset_j = 0.5 * (1 - decimal_i);
    C = [decimal_i - 1, decimal_j + offset_j - 0.5];
end

if INFO.is_south_pole
    TMP = N;
    N = S;
    S = TMP;
    N(1) = 4 * n - N(1);
    S(1) = 4 * n - S(1);
    W(1) = 4 * n - W(1);
    E(1) = 4 * n - E(1);
    C(1) = -C(1);
end


text:
classdef DWT2HD < dagnn.ElementWise

  properties (Transient)
      padding = 0
      wavename = 'haart'
      opts = {}
    numInputs
  end

  methods
    function outputs = forward(obj, inputs, params)
      obj.numInputs = numel(inputs) ;
      outputs{1} = vl_nndwt2(inputs{1}, [], ...
          'wavename', obj.wavename, 'padding', obj.padding, obj.opts{:}) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} =  vl_nndwt2(inputs{1},  derOutputs{1}, ...
          'wavename', obj.wavename, 'padding', obj.padding, obj.opts{:}) ;
      derParams = {0} ;
    end
    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes = {} ;
    end

    function rfs = getReceptiveFields(obj)
        rfs = [] ;
    end

    function obj = DWT2HD(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: MWCNN


text:
classdef cubochoricPlot < axisAnglePlot

  methods

    function oP = cubochoricPlot(varargin)
      % create a 3d plot of rotations in cubochoric coordinates

      oP = oP@axisAnglePlot(varargin{:});

     end

     function [x,y,z] = project(oP,ori,varargin)

      if ~check_option(varargin,'noBoundaryCheck')
        switch oP.fRMode
          case 'project2FundamentalRegion'
            ori = project2FundamentalRegion(ori);
          case 'restrict2FundamentalRegion'
            ori(~oP.oR.checkInside) = NaN;
        end
      end

      [x,y,z] = double(cubochoric(ori));

    end

    function ori = iproject(oP,x,y,z,varargin)
      ori = orientation.id;
    end

    function ori = makeGrid(oP,varargin)

      [ori,S2G,omega] = makeGrid@axisAnglePlot(oP,varargin{:});

      [oP.plotGrid.x, oP.plotGrid.y, oP.plotGrid.z] = ...
        double( S2G .* (3./4 * (omega - sin(omega))).^(1/3));

    end
  end
end


text:
function pass = test_max2(pref)
% Test the chebfun3/max2 command.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e7*pref.cheb3Prefs.chebfun3eps;

f = chebfun3(@(x,y,z) cos(x.*y.*z));
g = chebfun(@(x) 1 + 0*x);

h1 = max2(f);
h2 = max2(f, []);
h3 = max2(f, [], [1 2]);
h4 = max2(f, [], [2 1]);

h5 = max2(f, [], [1 3]);
h6 = max2(f, [], [3 1]);

h7 = max2(f, [], [2 3]);
h8 = max2(f, [], [3 2]);

pass(1) = norm(h1 - g) < tol;
pass(2) = norm(h2 - g) < tol;
pass(3) = norm(h3 - g) < tol;
pass(4) = norm(h4 - g) < tol;
pass(5) = norm(h5 - g) < tol;
pass(6) = norm(h6 - g) < tol;
pass(7) = norm(h7 - g) < tol;
pass(8) = norm(h8 - g) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function fim = fftconv2(im,f)
% function fim = fftconv2(im,f)
%
% Convolution using FFT.
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% March 2003

% wrap the filter around the origin and pad with zeros
padf = zeros(size(im));
r = floor(size(f,1)/2);
padf(1:r+1,1:r+1) = f(r+1:end,r+1:end);
padf(1:r,end-r+1:end) = f(r+2:end,1:r);
padf(end-r+1:end,1:r) = f(1:r,r+2:end);
padf(end-r+1:end,end-r+1:end) = f(1:r,1:r);

% magic
fftim = fft2(im);
fftf = fft2(padf);
fim = real(ifft2(fftim.*fftf));

author: zouchuhang

repo: LayoutNet


text:
function XY = rl_reward_firepit( type )
%RL_REWARD_FIREPITE
%
%   intput-----------------------------------------------------------------
%
%           o type : string, type of fire pit.
%
%   output-----------------------------------------------------------------
%
%           o XY   : (N x 2)
%
%

XY = [];

if strcmp(type,'central_pit')

    [X,Y] = meshgrid(linspace(3,7,20),linspace(3,7,20));
    XY    = [X(:),Y(:)];

else

end

end

author: epfl-lasa


text:
function data = seloverdim(data, seldim, sel, fb)

if nargin<4,
  fb = 1;
end

% get all XXXdimord fields
fn    = fieldnames(data);
selfn = find(~cellfun('isempty', strfind(fn, 'dimord')));
fn    = fn(selfn);
for k = 1:numel(fn)
  fndimord{k} = data.(fn{k});
end

% check which XXXdimord fields contain the seldim and keep only those
selx     = ~cellfun('isempty', strfind(fndimord, seldim));
fn       = fn(selx);
fndimord = fndimord(selx);

% extract the selected dimension as number
for k = 1:numel(fn)
  dimtok       = tokenize(fndimord{k}, '_');
  seldimnum{k} = find(strcmp(seldim, dimtok)); % the selected dimension as number
  if numel(seldimnum{k})<1 && strcmp(seldim, 'rpt'),
    %try rpttap
    seldimnum{k} = find(strcmp([seldim,'tap'], dimtok));
    tapflag      = 1;
  else
    tapflag      = 0;
  end
end

% if sum(~cellfun('isempty', seldimnum))<1
%   ft_error('the "%s" dimension is not present in the data', seldim)
% elseif any(cellfun(@numel, seldimnum)>1)
%   ft_error('cannot select over multiple dimensions at the same time')
% end

[reduceddim, fntmp] = dimlength(data);
selx       = find(ismember(fntmp, fn));
reduceddim = reduceddim(selx);
fntmp      = fntmp(selx);

% extract the fieldnames of the parameters of interest
if numel(fntmp)==1 && strcmp(fntmp{1}, 'dimord'),
  % data is not source data
  param      = selparam(data);
  reduceddim = repmat(reduceddim, [1 numel(param)]);
  seldimnum  = repmat(seldimnum,  [1 numel(param)]);
else
  for k = 1:numel(fntmp)
    param{1,k} = fntmp{k}(1:end-6);
  end
end

% make the subselection
for i = 1:numel(param)
  if numel(seldimnum{i})==1,
    if fb, fprintf('selection %s along dimension %d\n', param{i}, seldimnum{i}); end
  else
    if fb, fprintf('selection %s along dimensions %d and %d\n', param{i}, seldimnum{i}(1), seldimnum{i}(2)); end
  end

  reduceddim{i}(seldimnum{i}) = numel(sel);
  tmp       = data.(param{i});
  iscelltmp = iscell(tmp);
  if ~iscelltmp,
    %temporarily convert to cell
    tmp = {tmp};
  else
    %keep cells but reduce seldimnum
    seldimnum{i}  = seldimnum{i} - 1; %FIXME this only works if cell is 1D
    reduceddim{i} = [reduceddim{i}(2:end) 1];
  end

  if seldimnum{i}>0,
    % subselection from each of the cells
    for j = 1:numel(tmp)
      if ~isempty(tmp{j}) && numel(seldimnum{i})==1,
        switch seldimnum{i}
          case 1
            tmp{j} = tmp{j}(sel,:,:,:,:,:,:,:,:);
          case 2
            tmp{j} = tmp{j}(:,sel,:,:,:,:,:,:,:);
          case 3
            tmp{j} = tmp{j}(:,:,sel,:,:,:,:,:,:);
          case 4
            tmp{j} = tmp{j}(:,:,:,sel,:,:,:,:,:);
          case 5
            tmp{j} = tmp{j}(:,:,:,:,sel,:,:,:,:);
          case 6
            tmp{j} = tmp{j}(:,:,:,:,:,sel,:,:,:);
          case 7
            tmp{j} = tmp{j}(:,:,:,:,:,:,sel,:,:);
          case 8
            tmp{j} = tmp{j}(:,:,:,:,:,:,:,sel,:);
          case 9
            tmp{j} = tmp{j}(:,:,:,:,:,:,:,:,sel);
          otherwise
            ft_error('the number of dimensions is too high');
        end
        tmp{j} = reshape(tmp{j}, reduceddim{i});

      elseif ~isempty(tmp{j}) && numel(seldimnum{i})>1,
        if all(seldimnum{i}==[1 2])
          tmp{j} = tmp{j}(sel,sel,:,:,:,:,:,:,:,:);
        elseif all(seldimnum{i}==[2 3])
          tmp{j} = tmp{j}(:,sel,sel,:,:,:,:,:,:,:);
        else
          ft_error('selection of 2 dimensions simultaneously only works for [1 2] or [2 3] at present');
        end
        tmp{j} = reshape(tmp{j}, reduceddim{i});
      end
    end
  else
    % subselection of cells
    tmp = tmp(sel);
  end

  if ~iscelltmp,
    tmp = tmp{1};
  end
  data.(param{i}) = tmp;
end

switch seldim
  case 'rpt'
    if tapflag && isfield(data, 'cumtapcnt') && ~isempty(sel),
      sumtapcnt = cumsum([0;data.cumtapcnt(:)]);
      tapers    = zeros(1, sumtapcnt(end));
      for i=1:length(data.cumtapcnt)
        tapers(sumtapcnt(i)+1:sumtapcnt(i+1)) = i;
      end
      tmpsel    = [];
      tmpsel(1) = tapers(sel(1));
      for i=2:length(sel)
        if tapers(sel(i))~=tapers(sel(i-1))
	        tmpsel(end+1) = tapers(sel(i));
	      end
      end
    else
      tmpsel = sel;
    end
    if isfield(data, 'cumtapcnt'), data.cumtapcnt = data.cumtapcnt(tmpsel,:); end
    if isfield(data, 'cumsumcnt'), data.cumsumcnt = data.cumsumcnt(tmpsel,:); end

    % also try to adjust the trialinfo in the data
    if isfield(data, 'trialinfo')
      data.trialinfo = data.trialinfo(tmpsel, :);
    end

  case 'rpttap'
    % nothing to do
  case 'chan'
    data.label = data.label(sel);
  case 'freq'
    data.freq = data.freq(sel);
  case 'time'
    data.time = data.time(sel);
  case 'pos'
    data      = fixinside(data, 'logical');
    data.pos  = data.pos(sel,:);
    data.inside = data.inside(sel);
    data      = fixinside(data, 'index');
  otherwise
    ft_error('unknown dimension "%s"', seldim);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function kern = translateKernExpandParam(kern, params)

% TRANSLATEKERNEXPANDPARAM Create kernel structure from TRANSLATE kernel's parameters.
% FORMAT
% DESC returns a input space translation kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : translateKernParamInit, translateKernExtractParam,
% cmpndKernExpandParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

endVal = length(params)-kern.inputDimension;
kern = cmpndKernExpandParam(kern, params(1, 1:endVal));
startVal = endVal + 1;
kern.centre = params(1, startVal:end);


text:
function d2Vlims = opf_vlim_hess(x, lambda, mpc, idx, mpopt)
%OPF_VLIM_HESS  Evaluates Hessian of voltage magnitudes.
%   D2VLIMS = OPF_VLIM_HESS(X, LAMBDA, MPC, IDX, MPOPT)
%
%   Hessian evaluation function for voltage magnitudes.
%
%   Inputs:
%     X : optimization vector
%     LAMBDA : column vector of Lagrange multipliers on active and reactive
%              power balance constraints
%     MPC : MATPOWER case struct
%     IDX : index of buses whose voltage magnitudes should be fixed
%     MPOPT : MATPOWER options struct
%
%   Outputs:
%     D2VLIMS : Hessian of voltage magnitudes.
%
%   Example:
%       d2Vlims = opf_vlim_hess(x, lambda, mpc, idx, mpopt);
%
%   See also OPF_VLIM_FCN.

%   MATPOWER
%   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
%   by Baljinnyam Sereeter, Delft University of Technology
%   and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% unpack data
[Vr, Vi] = deal(x{:});

%% problem dimensions
nb = length(Vi);            %% number of buses
n = length(idx);            %% number of buses with voltage limits

%%----- evaluate Hessian of voltage limit constraints -----
nlam = length(lambda) / 2;
if nlam
    lam = lambda(nlam+1:2*nlam) - lambda(1:nlam);   %% lamVmax - lamVmin
else    %% keep dimensions of empty matrices/vectors compatible
    lam = zeros(0,1);
end

dlam = sparse(idx, idx, 2*lam, nb, nb);
zz = sparse(nb, nb);

%% construct Hessian
d2Vlims = [dlam zz; zz dlam];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function[fullpath]=findpath(dirname)
%FINDPATH  Returns the full pathname of a directory on the Matlab search path.
%
%   FINDPATH(DIRNAME) returns the full pathname of directory DIRNAME,
%   provided DIRNAME is on the Matlab search path.
%
%   For example, 'findpath jlab' on a Mac might return
%
%            '/Users/lilly/Desktop/DropBox/Matlab/jlab'.
%
%   FINDPATH returns the empty string if DIRNAME is not found.
%
%   Usage: findpath('dirname');
%          findpath dirname
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2015 J.M. Lilly --- type 'help jlab_license' for details

pathstr=path;
b=strfind(pathstr,':');
a=[1 b(1:end-1)+1];

fullpath=[];
for i=1:length(a)
    str=pathstr(a(i):b(i));
    index=strfind(str,[dirname ':']);
    if ~isempty(index)
        fullpath=str(1:end-1);
    end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function allFactors = BuildOCRNetwork (images, imageModel, pairwiseModel, tripletList)
% This function constructs the Markov network for one word by putting
% together the full list of different types of factors.
%
% Input:
%   images: A struct array of the images for each character in the word.
%   imageModel: The provided image model.
%   pairwiseModel: The provided pairwise factor model.
%   tripletList: The provided triplet factor list.
%
% Output:
%   allFactors: An array of factor structs that constitutes the full Markov
%     network.
%
% NOTE:
% The default implementation of each ConstructFactors function returns an
% empty list, so as you implement each of those in succession, more and
% more factors will be added to the resulting network.
%
% You can control which factors are added after you've implemented all of
% the factor functions in the following way:
%
% Pairwise factors: pass in an empty array (ie, []) for the pairwiseModel,
% and no pairwise factors will be added.
%
% Triplet factors: pass in an empty array for the tripletList, and no
% triplet factors will be added.
%
% Similarity factors: if you set imageModel.ignoreSimilarity = 1, then the
% similarity factors will not be added (nb: this is a bit of a kludge.
% Apologies.)
%
% Copyright (C) Daphne Koller, Stanford University, 2012

if (~exist('pairwiseModel', 'var'))
    pairwiseModel = [];
end
if (~exist('tripletList', 'var'))
    tripletList = [];
end

singletonFactors = ComputeSingletonFactors(images, imageModel);

if (~isempty(pairwiseModel))
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FOR STUDENTS:
    % Once you've tested out ComputeEqualPairwiseFactors, as discussed in
    % the assignment handout, you should comment out that line and
    % uncomment the following line to use the "real" pairwise factor
    % implementation (which you must provide.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%

   %pairwiseFactors = ComputeEqualPairwiseFactors(images, imageModel.K);
    pairwiseFactors = ComputePairwiseFactors(images, pairwiseModel, imageModel.K);
else
    pairwiseFactors = [];
end

if (~isempty(tripletList))
    tripletFactors = ComputeTripletFactors(images, tripletList, imageModel.K);
else
    tripletFactors = [];
end

% Your code here:

if (isfield(imageModel, 'ignoreSimilarity') && imageModel.ignoreSimilarity)
    simFactors = [];
else
    allSimFactors = ComputeAllSimilarityFactors(images, imageModel.K);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FOR STUDENTS:
    % In the provided implementation of ChooseTopSimilarityFactors, all
    % factors are returned, so no selection occurs. Once you try to run
    % inference using all factors (as described in the assignment handout),
    % you should implement the function so that the network contains only
    % the selected factors.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    simFactors = ChooseTopSimilarityFactors(allSimFactors, 2);
end

% This mess is necessary since Octave crashes if you try to vertcat onto a
% scalar, and it is possible that tripletFactors is scalar (for length 3
% words).
cellFactors{1} = singletonFactors;
cellFactors{2} = pairwiseFactors;
cellFactors{3} = tripletFactors;
cellFactors{4} = simFactors;
factorsPresent = ~[isempty(singletonFactors) isempty(pairwiseFactors) isempty(tripletFactors) isempty(simFactors)];

allFactors = vertcat(cellFactors{factorsPresent});
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function pdvfile( compleximagefile, outfile, varargin )
%PDVFILE Phase Derivative Value
%    pdvfile(compleximagefile, outfile, 'PropertyName',PropertyValue,...)
%    calculates the phase derivative value of compleximagefile using the
%    properties specified. This version of PDV does NOT require that the
%    complete data fit into memory.  It processes from any format handled
%    by OPEN_READER and outputs to files in SIO format.
%
%       Property name     Description
%       deltax            pixel shift (default = 0.25)
%       filtersize        size of smoothing filter (default = [5 5])
%       filtertype        type of filter, 'mean' (default) or 'median'
%       dim               dimension over which to calculate phase
%                            derivative (default = 1)
%       azlimits          min and max samples in azimuth over which to
%                            compute (default = query user with GUI).  Use
%                            'full' to specify entire range.  If nothing is
%                            selected, a MATLAB GUI will be used to select
%                            an area.
%       rnglimits         min and max lines in range over which to compute
%                            (default = query user with GUI).  Use 'full'
%                            to specify entire range.  If nothing is
%                            selected, a MATLAB GUI will be used to select
%                            an area.
%       framenumber      frame to process if a multi-image file
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Get image metadata
func_args=parsepdvinputs(varargin{:});
input_params=varargin;

% Read image metadata
fr=open_reader(compleximagefile);
if iscell(fr)
    metadata=fr{func_args.framenumber}.get_meta();
    for i=1:length(fr)
        fr{i}.close();
    end
else
    metadata=fr.get_meta();
    fr.close();
end

try
    % Normalize data if necessary.
    apply_normalization = ~is_normalized_sicd(metadata,func_args.dim);
    if apply_normalization
        tmp=tempname(); % Stored deskewed data in temporary file to process
        normalize_complex_file(compleximagefile,'output_filename',tmp,varargin{:});
        compleximagefile=tmp;
        input_params=[input_params {'azlimits','full','rnglimits','full'}];
    end
catch
    apply_normalization = 0;
end

% Actual PDV
func_hand=@(x) pdvmem(x, func_args);
process_by_lines( compleximagefile, outfile, func_hand,...
    'block_overlap', floor(func_args.filtersize.*[1 1]/2),...
    'function_name', 'PDV', input_params{:} );
if apply_normalization, delete(compleximagefile); end; % Remove temporary file

end

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

author: ngageoint


text:
function varargout = st_analyze(varargin)
% KJ_ANALYZE Application M-file for kj_analyze.fig
%    FIG = KJ_ANALYZE launch kj_analyze GUI.
%    KJ_ANALYZE('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.0 20-Aug-2003 09:47:18

if (nargin == 2) & ~ischar(varargin{1})  % LAUNCH GUI

  fig = openfig(mfilename);

  % Use system color scheme for figure:
  set(fig,'Color',get(groot,'defaultUicontrolBackgroundColor'));

  % Generate a structure of handles to pass to callbacks, and store it.
  handles = guihandles(fig);
  % Store the data
  handles.params = varargin{1};
  handles.hParentFigure = varargin{2};

  % Do a first plot with default values
  plotit(handles);

  guidata(fig, handles);

  if nargout > 0
    varargout{1} = fig;
  end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

  try
    [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
  catch
    disp(lasterr);
  end

end

% --------------------------------------------------------------------
function plotit(handles)

% Plot the frequency magnitude distribution
if get(handles.radEntireCatalog, 'Value') == 1
  mPlotCatalog = handles.params.mCatalog;
else
  mPlotCatalog = handles.mNodeCatalog;
end

axes(handles.axsFMD);
set(handles.axsFMD, 'NextPlot', 'replace');
set(handles.axsTimePlot, 'NextPlot', 'replace');
if get(handles.radEntireTime, 'Value') == 1
  % Plot the frequency magnitude distribution
  [hPlot, fBValue, fAValue, fStdDev, fMc] = plot_FMD(mPlotCatalog, 1, handles.axsFMD, 's', 'k', 1, handles.params.rOptions.nCalculateMC, handles.params.rOptions.fBinning);
  set(hPlot, 'MarkerSize', 10);
  sList = [cellstr('Entire time')];
  sList = [sList; cellstr(['     Mc: ' num2str(fMc)])];
  sList = [sList; cellstr(['b-value: ' num2str(fBValue)])];
  sList = [sList; cellstr(['a-value: ' num2str(fAValue)])];
  sList = [sList; cellstr(['std-dev: ' num2str(fStdDev)])];
  set(handles.lstValues, 'String', sList);
  % Plot the cumulative number
  [vDummy, vIndices] = sort(mPlotCatalog(:,3));
  mSortedCatalog = mPlotCatalog(vIndices(:,1),:) ;
  axes(handles.axsTimePlot)
  plot(mSortedCatalog(:,3),(1:length(mSortedCatalog(:,3))),'k-','era','xor');
else
  % Separate first catalog
  vSel_ = ((mPlotCatalog(:,3) >= handles.params.fStartFirstPeriod) & (mPlotCatalog(:,3) <= handles.params.fEndFirstPeriod));
  mFirst = mPlotCatalog(vSel_,:);

  % Separate second catalog
  vSel_ = ((mPlotCatalog(:,3) >= handles.params.fStartSecondPeriod) & (mPlotCatalog(:,3) <= handles.params.fEndSecondPeriod));
  mSecond = mPlotCatalog(vSel_,:);

  % Plot the frequency magnitude distribution for the first period
  [hPlot, fBValue1, fAValue, fStdDev, fMc] = plot_FMD(mFirst, 1, handles.axsFMD, 's', 'k', 1, handles.params.rOptions.nCalculateMC, handles.params.rOptions.fBinning);
  set(hPlot, 'MarkerSize', 3);
  sList = [cellstr('First period')];
  sList = [sList; cellstr(['     Mc: ' num2str(fMc)])];
  sList = [sList; cellstr(['b-value: ' num2str(fBValue1)])];
  sList = [sList; cellstr(['a-value: ' num2str(fAValue)])];
  sList = [sList; cellstr(['std-dev: ' num2str(fStdDev)])];
  sList = [sList; cellstr('')];
  % Plot the frequency magnitude distribution for the second period
  set(handles.axsFMD, 'NextPlot', 'add');
  [hPlot, fBValue2, fAValue, fStdDev, fMc] = plot_FMD(mSecond, 1, handles.axsFMD, '^', 'r', 1, handles.params.rOptions.nCalculateMC, handles.params.rOptions.fBinning);
  set(hPlot, 'MarkerSize', 3);
  sList = [sList; cellstr('Second period')];
  sList = [sList; cellstr(['     Mc: ' num2str(fMc)])];
  sList = [sList; cellstr(['b-value: ' num2str(fBValue2)])];
  sList = [sList; cellstr(['a-value: ' num2str(fAValue)])];
  sList = [sList; cellstr(['std-dev: ' num2str(fStdDev)])];
  sList = [sList; cellstr('')];
  sList = [sList; cellstr(['Delta-b: ' num2str(fBValue2-fBValue1)])];
  set(handles.lstValues, 'String', sList);

  % Normalize the two subcatalogs in time
  fMinTime = min(mFirst(:,3));
  mFirst(:,3) = mFirst(:,3) - fMinTime;
  fMinTime = min(mSecond(:,3));
  mSecond(:,3) = mSecond(:,3) - fMinTime;
  % Plot the cumulative number for the first period
  [vDummy, vIndices] = sort(mFirst(:,3));
  mSortedCatalog = mFirst(vIndices(:,1),:) ;
  axes(handles.axsTimePlot)
  plot(mSortedCatalog(:,3),(1:length(mSortedCatalog(:,3))),'k-','era','xor');
  % Plot the cumulative number for the second period
  [vDummy, vIndices] = sort(mSecond(:,3));
  mSortedCatalog = mSecond(vIndices(:,1),:) ;
  axes(handles.axsTimePlot)
  set(handles.axsTimePlot, 'NextPlot', 'add');
  plot(mSortedCatalog(:,3),(1:length(mSortedCatalog(:,3))),'r:','era','xor');
end

% --------------------------------------------------------------------
function varargout = radEntireCatalog_Callback(h, eventdata, handles, varargin)

if get(handles.radEntireCatalog, 'Value') == 1
  set(handles.radPickLocations, 'Value', 0);
  set(handles.btnPick, 'Enable', 'off');
else
  set(handles.radPickLocations, 'Value', 1);
  set(handles.btnPick, 'Enable', 'on');
end

% --------------------------------------------------------------------
function varargout = radPickLocations_Callback(h, eventdata, handles, varargin)

if get(handles.radPickLocations, 'Value') == 1
  set(handles.radEntireCatalog, 'Value', 0);
  set(handles.btnPick, 'Enable', 'on');
else
  set(handles.radEntireCatalog, 'Value', 1);
  set(handles.btnPick, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = radEntireTime_Callback(h, eventdata, handles, varargin)

if get(handles.radEntireTime, 'Value') == 1
  set(handles.radSplitPeriods, 'Value', 0);
else
  set(handles.radSplitPeriods, 'Value', 1);
end

% --------------------------------------------------------------------
function varargout = radSplitPeriods_Callback(h, eventdata, handles, varargin)

if get(handles.radSplitPeriods, 'Value') == 1
  set(handles.radEntireTime, 'Value', 0);
else
  set(handles.radEntireTime, 'Value', 1);
end

% --------------------------------------------------------------------
function varargout = btnRePlot_Callback(h, eventdata, handles, varargin)

plotit(handles);

% --------------------------------------------------------------------
function varargout = btnClose_Callback(h, eventdata, handles, varargin)

% Close the dialog
delete(handles.dlgAnalyze);

% --------------------------------------------------------------------
function varargout = btnPick_Callback(h, eventdata, handles, varargin)

% Get the axes handle of the plotwindow
axes(pf_result('GetAxesHandle', handles.hParentFigure, [], guidata(handles.hParentFigure)));

% Select a point in the plot window with the mouse
[fX, fY] = ginput(1);

pickit(handles, fX, fY);

% --------------------------------------------------------------------
function varargout = btnSelect_Callback(h, eventdata, handles, varargin)

fX = str2double(get(handles.txtXSelect, 'String'));
fY = str2double(get(handles.txtYSelect, 'String'));
pickit(handles, fX, fY);

% --------------------------------------------------------------------
function pickit(handles, fX, fY)

% Prepare the data
% If cross-section calculate the lenght along cross-section
if ~handles.params.bMap
  [nRow, nColumn] = size(handles.params.mCatalog);
  xsecx2 = handles.params.mCatalog(:,nColumn);  % length along x-section
  xsecy2 = handles.params.mCatalog(:,7);        % depth of hypocenters
end

% Calculate distance from center point and sort with distance
if handles.params.bMap
  vDistances = sqrt(((handles.params.mCatalog(:,1)-fX)*cos(pi/180*fY)*111).^2 + ((handles.params.mCatalog(:,2)-fY)*111).^2);
else
  vDistances = sqrt(((xsecx2 - fX)).^2 + ((xsecy2 + fY)).^2);
end

% Create the subcatalog
if handles.params.rOptions.nGriddingMode == 0
  % Use first nNumberEvents events
  [vTmp, vIndices] = sort(vDistances);
  mNodeCatalog = handles.params.mCatalog(vIndices,:);
  mNodeCatalog = mNodeCatalog(1:handles.params.nNumberEvents);
else
  % Use all events within fRadius
  vSelection = (vDistances <= handles.params.rOptions(1).fRadius);
  mNodeCatalog = handles.params.mCatalog(vSelection,:);
end

% Store the subcatalog for later use
handles.mNodeCatalog = mNodeCatalog;
guidata(handles.dlgAnalyze, handles);

% Export the node catalog
gui_export.mNodeCatalog = mNodeCatalog;
gui_export.fX = fX;
gui_export.fY = fY;
assignin('base', 'gui_export', gui_export);

% Replot the frequency magnitude distribution and cumulative number
plotit(handles);

% --------------------------------------------------------------------
function varargout = btnBWithTime_Callback(h, eventdata, handles, varargin)

axes(handles.axsBWithTime);
set(handles.axsBWithTime, 'NextPlot', 'replace');

%mCatalog = handles.mNodeCatalog;

fMcMin = str2double(get(handles.txtMcMin, 'String'));
fMcStep = str2double(get(handles.txtMcStep, 'String'));
fMcMax = str2double(get(handles.txtMcMax, 'String'));

nSteps = (fMcMax - fMcMin)/fMcStep + 1;
mColors = gui_Colormap_RedGrayBlue(nSteps);
nStep = 1;
for nMag = fMcMin:fMcStep:fMcMax
  vSel = handles.mNodeCatalog(:,6) >= nMag-0.0001;
  mCatalog = handles.mNodeCatalog(vSel,:);
% Defaults and inits
BV = [];
BV3 = [];
Nmin = 50;
nNumberEvents = str2double(get(handles.txtNumberEvents, 'String'));
nOverlap = str2double(get(handles.txtOverlap, 'String'));

for nCnt = 1:nNumberEvents/nOverlap:length(mCatalog) - nNumberEvents
    b = mCatalog(nCnt:nCnt + nNumberEvents,:);
    magco = min(b(:,6));
    l = b(:,6) >= magco-0.05;
    if length(b(l,:)) >= Nmin
        [uDummy bv stan,  uDummy] = bmemag(b(l,:));
    else
        bv = nan;
    end
    BV = [BV ; bv min(b(:,3)) ; bv max(b(:,3)) ; inf inf];
    BV3 = [BV3 ; bv mean(b(:,3)) stan ];
end

% Errorbar plot
hErrorbar = errorbar(BV3(:,2),BV3(:,1),BV3(:,3),BV3(:,3),'k');
set(hErrorbar(1), 'Color', [0.5 0.5 0.5], 'Linewidth', 1);
set(hErrorbar(2), 'Linewidth', 2, 'Color', mColors(nStep,:));
%
% % Timebar plot
% hTimebar = plot(BV(:,2),BV(:,1),'color',[0.5 0.5 0.5], 'Linewidth', 1);

% Values plot
%hValues = plot(BV3(:,2),BV3(:,1),'s');
%set(hValues,'LineWidth', 1, 'MarkerSize', 4,...
%    'MarkerFaceColor','w','MarkerEdgeColor','k','Marker','s');

  nStep = nStep + 1;
  set(handles.axsBWithTime, 'NextPlot', 'add');
end

author: CelsoReyes


text:
function d = Ddiv(A,B)
    d = sum(sum(A .* log (A./B) - A + B));
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef L1_Magic < handle
    % L1_Magic: using modified total variation of L1-magic package
    %
    % (c) Thomas Kuestner
    % ---------------------------------------------------------------------

    properties
%         oversampling % oversampling along space directions
        type % reconstruction type
        solver % define which optimization problem should be solved
        measPara % struct containing all measurement parameters (dim, dimension, ...)
%         dimension % 2D (x-y or x-y-t), 3D (x,y,z) or 4D (x,y,z,t) problem
%         dim % nPha x nFreq x nZ x nTime x nCha
        fullMask % mask of acquired data points
        trafo % struct of (tight frame) transformation in sparse basis
        epsilon % relaxation parameter
        lbtol % log barrier algorithm terminates when the duality gap <= lbtol
        mu % increase of barrier constant at each iteration
        cgtol % tolerance for Conjugate Gradients
        cgmaxiter % maximum number of iterations for Conjugate Gradients
        pdmaxiter % maximum number of iterations for primal dual algorithm
        pf % partial fourier parameters
        window % cut out window options
        kernelSize % kernel window dimension for convolution with calibration data
        currSlice % current slice, just needed for 2D with several slices
        currCha % current channel, just needed for 2D/3D with several channels
        lSense % apply SENSE mask for channel weighting
        kSpace % store kSpace in object
        espresso % ESPReSSo parameter
        dRecontime % recon time
    end

    methods
        function obj = L1_Magic(solver, measPara, epsilon, lbtol, mu, cgtol, cgmaxiter, pdmaxiter, espresso)
            obj.measPara = measPara;
            obj.epsilon = epsilon;
            obj.lbtol = lbtol;
            obj.mu = mu;
            obj.cgtol = cgtol;
            obj.cgmaxiter = cgmaxiter;
            obj.pdmaxiter = pdmaxiter;
            obj.solver = solver;
            obj.espresso = espresso;
            obj.type = 'L1_Magic';
        end
    end

end

author: thomaskuestner


text:
function AutoClustering(fbasename,elec,varargin)

% USAGE:
%     clu = AutoClustering(fbasename,elec)
%
% AutoClustering automtically cleans the clu file defined by fbasename and
% electrode numbers. The program will look for the files
% fbasename.fet/res/spk.elec in the current folder.
%
% INPUT:
%     fbasename: char array
%     elec: a vector of electrode numbers
%
% optional:
%     AutoClustering(fbasename,elec,dim)
%     where dim is the number of channels in electro group (if not
%     defined, will read the first line of the fet file
%
% AutoClustering is meant to clean the output of KlustaKwik. The first
% thing it does is to separate electrical artifacts and MUA from putative
% isolated units. To do so, it sorts out units which have no clear
% refractory period (based on Hill, Mehta and Kleinfeld, J Neurosci.,
% 2012). Threshold can be set in the parameter section of this file
% ("Rogue spike threshold"). Then, it separates electrical
% artifats from MUA based on the assumption that electrical artifacts are
% highly correlated on thel different channels: the average waveform of at
% least one channel has to be different from the across-channel average
% waveform by a certrain amount of total variance (can be set in the
% parameter section, "Deviation from average spike threshold")
%
% Once the program has determined which of the clusters are putative
% isolated units, it tries to merge them based on waveform similarity
% (mahalanobis distance) and quality of the refractory period in the new
% merged cluster (or "Inter Common Spike Interval" from MS Fee et al. J
% Neurosci. Meth., 1996)
%
% Adrien Peyrache, 2012
% David Tingley, 2017 updated for Intan system and new functionality

% The default behavior of this program is to take the output from KlustaKwik
% or the newer klusta-3.0 (.kwik) algorithms, organize this data for easier
% manual cluster cutting, and save to the klusta-3.0 (.kwik) format

if ~isempty(dir('*kwik'))
 if nargin < 1
    basepath = pwd;
    name = dir('*alg*');
    s = split(name.name,'.');
    fbasename = s{1};
    elec = str2num(s{end});
 end

dbstop if error
% Parameters
% Number recording sites
if ~isempty(varargin)
    dim = varargin{1};
    dim = dim(:);
    if any(double(int16(dim))~=dim)
        error('Number of dimensions must be an integer')
    end
    if size(dim,1) ~= length(elec) && length(dim) ~=1
        error('Number of dimensions must be a vector of the same length as electrode vecotr or a single value')
    end
    if length(dim) == 1
        dim = dim*ones(length(elec),1);
    end
else
    dim = zeros(length(elec),1);
end

% an ugly list of heuristics that should eventually be cleaned up...

samplingRate = 20000;
% Load spike waveforms? Used for detection of electrical artifacts
loadspk = 0;
% Refractory period in msec
tR = 1.5./1000;
% Censored period in msec (specific to the spike detection process)
tC = 0.01./1000;
% Rogue spike threshold (for MUA); value between 0 an 1
rogThres = 1.5;
% Relative deviation (from 0 to 1) from average spike threshold (for electrical artifacts)
% =1000 => bypass it
devMinThres = .3;
devMaxThres = 10;

% other arbitrary thresholds..
powThresh = .85;
isoMinTresh = 1;
lratioMinThresh = 10e-5;
isoMaxTresh = 2000;
lratioMaxThresh = 1000;

% Do Merging?
doMerge = 1;
% overwrite clu file
rewriteclu= 1;
% Write a log file?
WriteLogFile = 1;

if WriteLogFile
    tic;
    log = [];
end

elec = elec(:)';
if length(elec)>1
    for eIx=1:length(elec)
        AutoClustering(fbasename,elec(eIx),dim(eIx));
    end
else
    % Load fet, clu, res, and waveforms
    fprintf('Sorting electrode %i of %s\n',elec,fbasename)
    if 1 %~exist([fbasename '_sh' num2str(elec) '.res.' num2str(elec)]) & ~exist([fbasename '.res.' num2str(elec)])
        if exist([fbasename '_sh' num2str(elec) '.kwik']) > 0
            tkwik = fullfile(pwd,[fbasename '_sh' num2str(elec) '.kwik']);
            cd ..; basepath = pwd; cd(num2str(elec));
        else
            tkwik = fullfile(pwd,num2str(elec),[fbasename '_sh' num2str(elec) '.kwik']);
            basepath = pwd;
        end
        kwikinfo = h5info(tkwik,['/channel_groups/' num2str(elec) '/clusters/original']);
        if ~loadspk
            [fet clu res ~] = ConvertKlusta2Matlab(elec,basepath,fbasename,0,0,1);
        elseif loadspk
           [fet clu res wav] = ConvertKlusta2Matlab(elec,basepath,fbasename,1,0,1);
        end
        clu_orig = clu;
        clu = double(clu);
        cluster_names = unique(clu);
    else
        error('check the below code for compatibility...')
        fet = LoadFeatures(fbasename,elec,dim);
        clu = load([fbasename '.clu.' num2str(elec)]);
        clu = clu(2:end);
        res = load([fbasename '.res.' num2str(elec)])/20;
        if loadspk
            wav = LoadSpikeWaveforms([fbasename '.spk.' num2str(elec)],dim,32);
        end
    end

    %Power of each waveform - UNUSED so far
%     pw = sqrt(sum(sum(wav.^2,2),3));
    % Percentile of the power
%     p = prctile(pw,0.99);
    nFeats = size(fet,2);  % to calculate accurate stats, we should have more rows than features for all cells

    %  devFromMean quantifies how much the spike are different on
    %  the different channels (what is the maximal distance from the averaged
    %  spike)

    for ii=1:length(cluster_names)
        rg = double(res(clu==cluster_names(ii)))./samplingRate;
        if loadspk
            m = squeeze(mean(wav(clu==cluster_names(ii),:,:)));
            y = pdist(m,'euclidean')./max(sqrt(sum(m.^2,2)));
            devFromMean(ii) = max(y);
        end
        fractRogue(ii) = FractionRogueSpk(rg,tR,tC);
%         meanPw(ii) = nanmean(pw(clu==cluster_names(ii)))/nanmean(pw);
    end
    if loadspk
        ff = find([devFromMean < devMinThres |...
             devFromMean > devMaxThres]);  % enforces variability across channels (muscle)
    else
        ff = [];
    end

    fff = find([fractRogue > rogThres]);  % enforces refractory period

    for ii=1:length(cluster_names)
        if length(find(clu==cluster_names(ii))) > nFeats
            [L(ii) LRatio(ii)] = L_Ratio(fet,find(clu==cluster_names(ii)));
            iso(ii) = IsolationDistance(fet,find(clu==cluster_names(ii)));
        else
            L(ii)=nan;LRatio(ii)=nan;iso(ii)=nan; % not enough spikes to quantify anything...
        end
    end
    % Here let's toss out clusters that are clearly muscle/electrical
    % we have: fractRogue, devFromMean, LRatio, and iso to use
    f = find([LRatio<lratioMinThresh |...
            LRatio > lratioMaxThresh &...
            iso < isoMinTresh |...
            iso > isoMaxTresh]);          % finds large artefacts
    noiseIx = find(ismember(clu,cluster_names(unique([f,ff,fff]))));

    % takes all muscle/electrical artefact and merges to a single
    % garbage cluster
    badCluIdx = unique(clu(noiseIx));
    for ii=1:length(badCluIdx)
        log = [log sprintf('%d -> %d; Looks like muscle or electrical artifact\n',badCluIdx(ii),0)];
    end
    clu(noiseIx) = 0;
    cluster_names = unique(clu);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Display intermediary results
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if 0
        for i=1:length(cluster_names)
            subplot(2,2,1)
            m = squeeze(mean(wav(clu==cluster_names(i),:,:)));
            plot(m'); title(['Deviation from mean wf: ' num2str(devFromMean(i))])
            subplot(2,2,2)
            rg = double(res(clu==cluster_names(i)))./samplingRate;
            dt = diff(rg);
            hist(dt(dt<.1),1000); title(['Fraction Rogue: ' num2str(fractRogue(i))]);
            subplot(2,2,3)
            title(['L-Ratio: ' num2str(LRatio(i)) ' IsoDist: ' num2str(iso(i))])
            subplot(2,2,4)
            if ismember(cluster_names(i),unique(clu(noiseIx)))
                title('tossing');
            else
                title('keeping');
            end
%             title(num2str(d(:,i)'))
            pause
        end
    end

    %reorder clusters here
    [clu log] = renumberclu(clu,log);

    % Here we select only clusters that correspond to putative units and that
    % have at least 20 spikes (otherwise errormatrix calculation fails)
    h = hist(clu,length(unique(clu)))';
    goodCluIx = ismember(clu,find(cluster_names < 1000 & h>nFeats));
    goodCluIx(noiseIx) = 0;

    if length(unique(clu(goodCluIx)))>1  % if we have more than one putative cluster...
        newclu = clu(goodCluIx);
        newres = res(goodCluIx);
        newfet = fet(goodCluIx,:);
        if doMerge  % merge similar clusters which are neither noise nor MUA
            try
                [newclu mergehistory] = mergeclu_slow(newclu,newres,newfet,tR,tC,rogThres,samplingRate);
                % log changes
                log = [log sprintf('merge_slow.m was run\n')];
                for ii=1:size(mergehistory,1)
                    log = [log sprintf('%d + %d -> %d\n',mergehistory(ii,1),mergehistory(ii,2),mergehistory(ii,3))];
                end
            catch
                warning(['merging failed: ' lasterr])
            end
        else   % if we're not going to try merging, the resort clusters by similarity
            em = errormatrix(fet(goodCluIx,:),newclu);
            ems = max(em,em');
            y = squareform((1-ems)-diag(diag(1-ems)),'tovector');
            Z = linkage(y);
            [T,perm] = dendrogram_struct(Z,0);
            cluIx = unique(newclu);
            for ii=1:length(cluIx)
                newclu = updateclu(newclu,cluIx(perm(ii)),max(cluIx)+ii+1);
            end
        end
        clu(goodCluIx) = newclu; % re-insert any merges or resorting of cluster ID's
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Display final results (requires function CrossCorr, not in the
    % toolbox)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if 0
    cluIx = unique(newclu(newclu>1));
    for x=1:length(cluIx)
        rgx = res(newclu==cluIx(x));
        [ax,b] =  CrossCorr(rgx*10,rgx*10,1,60);ax(b==0)=0;
        for y=x+1:length(cluIx)
            rgy = res(newclu==cluIx(y));
            [ay,b] =  CrossCorr(rgy*10,rgy*10,1,60);ay(b==0)=0;
            [fxy fcorr] = crossrefract(rgx,rgy);
            [h,b] = CrossCorr(rgx*10,rgy*10,1,60);ac(b==0)=0;
            figure(1),clf
            subplot(1,3,1)
                bar(b,ax,1)
             subplot(1,3,2)
                bar(b,h,1)
                title([fxy fcorr])
             subplot(1,3,3)
                bar(b,ay,1)
            pause
        end
    end
    end

     % reorder clusters here
    [clu log] = renumberclu(clu,log);
    disp(log)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Write new clu file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if rewriteclu
        kwikinfo = h5info(tkwik,['/channel_groups/' num2str(elec) '/clusters/main']);
        kwikinfo_original = h5info(tkwik,['/channel_groups/' num2str(elec) '/clusters/original']);

        if length(cluster_names) > length(kwikinfo_original.Groups)
%             kwikinfo = h5info(tkwik,['/channel_groups/' num2str(elec) '/clusters/main']);
            error('we added to the number of clusters?')
        end
        % rewrite cluster ID's
        h5write(tkwik,['/channel_groups/' num2str(elec) '/spikes/clusters/main'],uint32(clu));

%         I need to change autoclustering to selectively delete cluster groups that no longer have units....
        % delete all old cluster gropus from /main
%         fid = H5F.open(kwikinfo.Filename,'H5F_ACC_RDWR','H5P_DEFAULT');
%         for gg = 1:length(kwikinfo.Groups)
%             H5L.delete(fid,[kwikinfo.Groups(gg).Name],'H5P_DEFAULT')
%         end
%         H5F.close(fid)
        % copy cluster groups from /original to /main for new clusters
        fid = H5F.open(tkwik,'H5F_ACC_RDWR','H5P_DEFAULT');
        for i =1:length(cluster_names) % rewrite cluster groups
            try  % this may fail if something didn't get deleted... that's ok
            H5L.copy(fid,kwikinfo_original.Groups(i).Name,...
                            fid,[kwikinfo.Name '/' num2str(cluster_names(i))],...
                            'H5P_DEFAULT','H5P_DEFAULT')
            catch
                disp([ num2str(cluster_names(i)) ' cluster group already exists'])
            end
            if doMerge
                h5writeatt(tkwik,['/channel_groups/' num2str(elec) '/clusters/main/'...
                    num2str(cluster_names(i))],'cluster_group',3) % a human has not looked at these results, everything should be 'unsorted' (3 in .kwik format)
            end
        end
        H5F.close(fid)

        if exist([num2str(elec) '/nohup.out'])  && exist([num2str(elec) '/' fbasename '_sh' num2str(elec) '.klg.' num2str(elec)])
            if exist([fbasename '_sh' num2str(elec) '.kwik']) > 0
                fileID = fopen('nohup.out','a');
                klgID = fopen([fbasename '_sh' num2str(elec) '.klg.' num2str(elec)],'a');
            else
                fileID = fopen([num2str(elec) '/nohup.out'],'a');
                klgID = fopen([num2str(elec)  '/' fbasename '_sh' num2str(elec) '.klg.' num2str(elec)],'a');
            end
            fmt = 'this elec has been autoclustered';
            fprintf(fileID,fmt);
            fprintf(klgID,fmt);
            fclose(klgID);
            fclose(fileID); % write to both nohup and .klg. log files
        else
            klgID = fopen([fbasename '_sh' num2str(elec) '.klg.' num2str(elec)],'a');
            fmt = 'this elec has been autoclustered';
            fprintf(klgID,fmt);
            fclose(klgID);
            warning('could not find nohup.out log file')
        end
    end

    if WriteLogFile
        % Create (of overwrite) a log file
        fid = fopen([fbasename '.alg.' num2str(elec)],'w');
        time_tot = toc;
        fprintf(fid,[log 'That took %f seconds\n'],time_tot);
    end
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function indx = r8col_sort_heap_index_a ( m, n, a )

%*****************************************************************************80
%
%% R8COL_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R8COL.
%
%  Discussion:
%
%    The sorting is not actually carried out.  Rather an index array is
%    created which defines the sorting.  This array may be used to sort
%    or index the array, or to sort or index related arrays keyed on the
%    original array.
%
%    A(*,J1) < A(*,J2) if the first nonzero entry of A(*,J1)-A(*,J2) is negative.
%
%    Once the index array is computed, the sorting can be carried out
%    "implicitly:
%
%      A(*,INDX(1:N)) is sorted,
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in each column of A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the array.
%
%    Output, integer INDX(N), the sort index.  The I-th element of the sorted
%    array is column INDX(I).
%
  if ( n < 1 )
    indx = [];
    return
  end

  indx(1:n) = 1 : n;

  if ( n == 1 )
    return
  end

  l = floor ( n / 2 ) + 1;
  ir = n;

  while ( 1 )

    if ( 1 < l )

      l = l - 1;
      indxt = indx(l);
      column(1:m) = a(1:m,indxt);

    else

      indxt = indx(ir);
      column(1:m) = a(1:m,indxt);
      indx(ir) = indx(1);
      ir = ir - 1;

      if ( ir == 1 )
        indx(1) = indxt;
        break
      end

    end

    i = l;
    j = l + l;

    while ( j <= ir )

      if ( j < ir )

        if ( r8vec_compare ( m, a(1:m,indx(j)), a(1:m,indx(j+1)) ) < 0 )
          j = j + 1;
        end

      end

      if ( r8vec_compare ( m, column, a(1:m,indx(j)) ) < 0 )
        indx(i) = indx(j);
        i = j;
        j = j + j;
      else
        j = ir + 1;
      end

    end

    indx(i) = indxt;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [model, B, elapse] = USPLH_learn(A, maxbits)
%   This is a wrapper function of USPLH learning.
%
%	Usage:
%	[model, B,elapse] = USPLH_learn(A, maxbits)
%
%	      A: Rows of vectors of data points. Each row is sample point
%   maxbits: Code length
%
%     model: Used for encoding a test sample point.
%	      B: The binary code of the input data A. Each row is sample point
%    elapse: The coding time (training time).
%
%
%
%   version 2.0 --Nov/2016
%   version 1.0 --Jan/2013
%
%   Written by  Yue Lin (linyue29@gmail.com)
%               Deng Cai (dengcai AT gmail DOT com)
%

tmpT = tic;

model = trainUSPLH(A, struct('nbits',maxbits,'eta',0.125));
B = USPLH_compress(A, model);

elapse = toc(tmpT);

end

author: ZJULearning

repo: MatlabFunc


text:
function zvalues = interactivedendrogramvalues(fig)
% zvalues = interactivedendrogramvalues(fig)
% from http://www.mathworks.com/help/techdoc/ref/datacursormode.html#bsawkea-7
dcm_obj = datacursormode(fig);
stop = 0;
inp = input('How many nodes?\n');
for ii=1:inp

    set(dcm_obj,'DisplayStyle','datatip',...
        'SnapToDataVertex','off','Enable','on')

    disp('Click on a node, then press Return.')
    figure(fig)
    pause                            % Wait while the user does this.
    c_info = getCursorInfo(dcm_obj);
    zvalues(ii)=c_info.Position(2);
    fprintf('Zvalue:%g\n', zvalues(ii))
end

% while ~stop
%     set(dcm_obj,'DisplayStyle','datatip',...
%         'SnapToDataVertex','off','Enable','on')
%
%     disp('Click on a node, then press Return.')
%     pause                            % Wait while the user does this.
%
%     c_info = getCursorInfo(dcm_obj);
%     zvalues(ii)=c_info.Position(2);
%     ii=ii+1;
%
%     inp = input('More? Yes = [anything but n]/ No = [n]\n','s');
%
%     if strcmp(inp,'n')
%         stop=1;
%     end
% end


text:
clear;  clc;
%P2.16
%y(n)-0.5y(n-1)+0.25y(n-2)=x(n)+2x(n-x)+x(n-3)
a=[1,-0.5,0.25];b=[1,2,1];
h=impseq(0,-10,100);n=[-10:100];
y=filter(b,a,h);
subplot(2,1,1);
stem(n,y,'k');title('Impulse Response');xlabel('n');ylabel('h(n)')
%We determine S|h(n)|
sum(abs(y))
magz = abs(roots(a))
%ans = 6.5714
%Which implies that the system is stable.

%If x(n)=[5+cos(0.2pn)+4sin(0.6pn)]u(n),
%the impulse response will show no stability
a=[1,-0.5,0.25];b=[1,2,1];
n=[-10:200];
x=(5+3*cos(0.2*pi*n)+4*sin(0.6*pi*n)).*stepseq(0 ,-10,200);
y=filter(b,a,x);
subplot(2,1,2);
stem(n,y,'k');title('Response to x(n)');xlabel('n');ylabel('y(n)')
sum(abs(y))
magz = abs(roots(a))
%ans =  5.3389e+003
% We observe that the response is not bounded, since the output is not close to zero as n approaches infinity.

author: Sable


text:
function prob_test012 ( )

%*****************************************************************************80
%
%% PROB_TEST012 tests BETA_CDF, BETA_CDF_INV, BETA_PDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROB_TEST012\n' );
  fprintf ( 1, '  For the Beta PDF:\n' );
  fprintf ( 1, '  BETA_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  BETA_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  BETA_PDF evaluates the PDF;\n' );

  a = 12.0;
  b = 12.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );

  check = beta_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PROB_TEST012 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '        A               B               X               ' )
  fprintf ( 1, 'PDF             CDF             CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = beta_sample ( a, b, seed );

    pdf = beta_pdf ( x, a, b );

    cdf = beta_cdf ( x, a, b );

    x2 = beta_cdf_inv ( cdf, a, b );

    fprintf ( 1, '%14f  %14f  %14f  %14f  %14f  %14f\n', a, b, x, pdf, cdf, x2 );

  end

  return
end


text:
classdef ScalarGaussVariablePrinter < ResultsPrinter

    properties (Access = protected)
        simulationStr
        hasGaussData = true;
    end

    properties (Access = private)
        fieldName
    end

    methods (Access = public)

        function obj = ScalarGaussVariablePrinter(d)
            obj.init(d);
            obj.fieldName = d.fieldName;
            obj.simulationStr = [d.fieldName,'Gauss'];
        end

        function printResults(obj,iter,fileID)
            field(:,1,:) = squeeze(obj.fields)';
            dS = obj.createScalarGaussDataBase(iter,fileID,field, obj.fieldName,'OnGaussPoints');
            ScalarGaussPrinter(dS);
        end

    end

    methods (Access = protected)

        function storeFieldsToPrint(obj,d)
            obj.fields = d.fields;
        end

        function createHeadPrinter(obj,d,dh)
            obj.headPrinter = GaussHeadPrinter(d,dh);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = imageProc(varargin)

% IMAGEPROC M-file for imageProc.fig
%      IMAGEPROC, by itself, creates a new IMAGEPROC or raises the existing
%      singleton*.
%
%      H = IMAGEPROC returns the handle to a new IMAGEPROC or the handle to
%      the existing singleton*.
%
%      IMAGEPROC('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in IMAGEPROC.M with the given input arguments.
%
%      IMAGEPROC('Property','Value',...) creates a new IMAGEPROC or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before imageProc_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to imageProc_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Copyright 2002-2003 The MathWorks, Inc.

% Edit the above text to modify the response to help imageProc

% Last Modified by GUIDE v2.5 15-Feb-2011 10:30:37

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @imageProc_OpeningFcn, ...
                   'gui_OutputFcn',  @imageProc_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before imageProc is made visible.
function imageProc_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to imageProc (see VARARGIN)

% Choose default command line output for imageProc
handles.fileLoaded = 0;
handles.fileLoaded2 = 0;
set(handles.axes1,'Visible','off');
set(handles.axes2,'Visible','off');
set(handles.axesHist1,'Visible','off');
set(handles.axesHist2,'Visible','off');
set(handles.editPath, 'Visible', 'off');
set(handles.editSize, 'Visible', 'off');
set(handles.editComment, 'Visible', 'off');
set(handles.textHist1, 'Visible', 'off');
set(handles.textHist2, 'Visible', 'off');
set(handles.sliderBright, 'Enable', 'off');
set(handles.sliderContrast, 'Enable', 'off');
set(handles.sliderRotate, 'Enable', 'off');
set(handles.editBright,'String', sprintf('%10s:%4.0f%%', 'Brightness', 100*get(handles.sliderBright,'Value')));
set(handles.editContrast,'String', sprintf('%10s:%4.0f%%', 'Contrast', 100*get(handles.sliderContrast,'Value')));
set(handles.editRotate,'String', sprintf('%10s:%4.0f', 'Rotate', get(handles.sliderRotate,'Value')));
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes imageProc wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = imageProc_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in LoadButton.
function LoadButton_Callback(hObject, eventdata, handles)
% hObject    handle to LoadButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[FileName,PathName] = uigetfile({'*.*'},'Load Image File');

if (FileName==0) % cancel pressed
    return;
end

handles.fullPath = [PathName FileName];
[a, b, Ext] = fileparts(FileName);
availableExt = {'.bmp','.jpg','.jpeg','.tiff','.png','.gif'};
FOUND = 0;
for (i=1:length(availableExt))
    if (strcmpi(Ext, availableExt{i}))
        FOUND=1;
        break;
    end
end

if (FOUND==0)
    h = msgbox('File type not supported!','Error','error');
    return;
end

set(handles.sliderRotate, 'Enable', 'on');
set(handles.sliderBright, 'Enable', 'on');
set(handles.sliderContrast, 'Enable', 'on');
set(handles.editPath, 'Visible', 'on');
set(handles.editSize, 'Visible', 'on');
set(handles.editComment, 'Visible', 'on');

info = imfinfo(handles.fullPath);
if (~isempty(info.Comment))
    % save current image comment (to be used later in image save)
    handles.currentImageComment = info.Comment{1};
else
    handles.currentImageComment = '';
end

set(handles.editSize, 'String', sprintf('SIZE (W x H) : %d x %d', info.Width, info.Height));
set(handles.editComment, 'String', sprintf('COMMENT: %s', handles.currentImageComment));
set(handles.editPath', 'String', handles.fullPath);

RGB = imread(handles.fullPath);

handles.RGB = RGB;
handles.RGB2 = RGB;
handles.fileLoaded = 1;
handles.fileLoaded2 = 0;

set(handles.axes1,'Visible','off'); set(handles.axes2,'Visible','off');
set(handles.axesHist1,'Visible','off'); set(handles.axesHist2,'Visible','off');
set(handles.textHist1, 'Visible', 'off');
axes(handles.axesHist2); cla;
set(handles.textHist2, 'Visible', 'off');

axes(handles.axes1); cla; imshow(RGB);
axes(handles.axes2); cla;

handles = updateHistograms(handles);

guidata(hObject, handles);

% --- Executes on button press in CopyButton.
function CopyButton_Callback(hObject, eventdata, handles)
% hObject    handle to CopyButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded==1)
    handles.RGB2 = handles.RGB;
    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);
    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

% --- Executes on button press in MedianButton.
function MedianButton_Callback(hObject, eventdata, handles)
% hObject    handle to MedianButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded==1)

    [M,N,ttt] = size(handles.RGB);

    RUN = 1;

    while (RUN==1)

        prompt = {'Enter Median Row Factor (0-5%):','Enter Median Column Factor (0-5%):'};
        dlg_title = 'Enter Median Parameters:';
        num_lines = 1;
        def = {'2','2'};
        answer = inputdlg(prompt,dlg_title,num_lines,def);
        if (isempty(answer))
            return;
        end

        M1 = str2num(answer{1})/100;
        M2 = str2num(answer{2})/100;

        if ((str2num(answer{1})>=0) & (str2num(answer{1})<=5)) & ((str2num(answer{2})>=0) & (str2num(answer{2})<=5))
            RUN = 0;
        end
    end

    M1 = round(M1 * M);
    M2 = round(M2 * N);

    w = waitbar(0, 'Median filtering ... Please wait ...');
    handles.RGB2(:,:,1) = medfilt2(handles.RGB(:,:,1),[M1 M2]);
    waitbar(1/3, w);
    handles.RGB2(:,:,2) = medfilt2(handles.RGB(:,:,2),[M1 M2]);
    waitbar(2/3, w);
    handles.RGB2(:,:,3) = medfilt2(handles.RGB(:,:,3),[M1 M2]);
    close(w);
    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);

    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

% --- Executes on button press in SharpButton.
function SharpButton_Callback(hObject, eventdata, handles)
% hObject    handle to SharpButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded==1)
    H = fspecial('unsharp');
    handles.RGB2(:,:,1) = imfilter(handles.RGB(:,:,1),H,'replicate');
    handles.RGB2(:,:,2) = imfilter(handles.RGB(:,:,2),H,'replicate');
    handles.RGB2(:,:,3) = imfilter(handles.RGB(:,:,3),H,'replicate');
    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);
    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

% --- Executes on button press in MotionButton.
function MotionButton_Callback(hObject, eventdata, handles)
% hObject    handle to MotionButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded==1)

    [M,N,ttt] = size(handles.RGB);
    D = max(M,N);

    RUN = 1;

    while (RUN==1)

        prompt = {'Enter Motion Length (0-15%):','Enter Motion Angle (0-360):'};
        dlg_title = 'Enter Motion Parameters:';
        num_lines = 1;
        def = {'1','0'};
        answer = inputdlg(prompt,dlg_title,num_lines,def);
        if (isempty(answer))
            return;
        end

        M1 = str2num(answer{1})/100;
        M2 = str2num(answer{2});

        if ((str2num(answer{1})>=0) & (str2num(answer{1})<=15))
            RUN = 0;
        end
    end

    H = fspecial('motion',D * M1,M2);
    w = waitbar(0, 'Motion filtering ... Please wait ...');
    handles.RGB2(:,:,1) = imfilter(handles.RGB(:,:,1),H,'replicate');
    waitbar(1/3, w);
    handles.RGB2(:,:,2) = imfilter(handles.RGB(:,:,2),H,'replicate');
    waitbar(2/3, w);
    handles.RGB2(:,:,3) = imfilter(handles.RGB(:,:,3),H,'replicate');
    close(w);
    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);
    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

% --- Executes on button press in GrayButton.
function GrayButton_Callback(hObject, eventdata, handles)
% hObject    handle to GrayButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded==1)
    Gray = rgb2gray(handles.RGB);
    handles.RGB2(:,:,1) = Gray;
    handles.RGB2(:,:,2) = Gray;
    handles.RGB2(:,:,3) = Gray;
    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);
    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

% --- Executes on button press in SaveButton.
function SaveButton_Callback(hObject, eventdata, handles)
% hObject    handle to SaveButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded2==1)
    [file,path] = uiputfile('*.jpg','Save Secondary Image As');
    imwrite(handles.RGB2,[path file],'jpg');
else
    h = msgbox('No secondary file has been loaded!','Save Error','error');
end

% --- Executes on button press in ColorsButton.
function ColorsButton_Callback(hObject, eventdata, handles)
% hObject    handle to ColorsButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded==1)
    RUN = 1;

    while (RUN==1)

        prompt = {'Enter threshold for RED (0-255):','Enter threshold for GREEN (0-255):','Enter threshold for BLUE (0-255):'};
        dlg_title = 'RGB Thresholds:';
        num_lines = 1;
        def = {'30','30','30'};
        answer = inputdlg(prompt,dlg_title,num_lines,def);
        if (isempty(answer))
            return;
        end
        T1 = str2num(answer{1});
        T2 = str2num(answer{2});
        T3 = str2num(answer{3});

        if ((T1>=0) & (T1<=256)) & ((T2>=0) && (T2<=256)) & ((T3>=0) && (T3<=256))
            RUN = 0;
        end
    end

    handles.RGB2 = filterColors(handles.RGB, T1, T2, T3, 5);
    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);

    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

% --- Executes on button press in ColorButton2.
function ColorButton2_Callback(hObject, eventdata, handles)
% hObject    handle to ColorButton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded==1)
    RUN = 1;

    while (RUN==1)

        prompt = {'Enter weight for RED (0-200%):','Enter weight for GREEN (0-200%):','Enter weight for BLUE (0-200%):'};
        dlg_title = 'Enter Color Weight Parameters:';
        num_lines = 1;
        def = {'100','100','100'};
        answer = inputdlg(prompt,dlg_title,num_lines,def);
        if (isempty(answer))
            return;
        end
        W1 = str2num(answer{1});
        W2 = str2num(answer{2});
        W3 = str2num(answer{3});

        if ((W1>=0) & (W1<=200)) & ((W2>=0) && (W2<=200)) & ((W3>=0) && (W3<=200))
            RUN = 0;
        end
    end

    handles.RGB2 = handles.RGB;
    R = double(handles.RGB2(:,:,1));
    G = double(handles.RGB2(:,:,2));
    B = double(handles.RGB2(:,:,3));

    R = R * W1 ./ 100;
    G = G * W2 ./ 100;
    B = B * W3 ./ 100;

    R(find(R>256)) = 256;
    G(find(G>256)) = 256;
    B(find(B>256)) = 256;

    handles.RGB2(:,:,1) = R;
    handles.RGB2(:,:,2) = G;
    handles.RGB2(:,:,3) = B;

    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);
    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

% --- Executes on button press in InvColorButton.
function InvColorButton_Callback(hObject, eventdata, handles)
% hObject    handle to InvColorButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if (handles.fileLoaded==1)

    R = double(handles.RGB(:,:,1));
    G = double(handles.RGB(:,:,2));
    B = double(handles.RGB(:,:,3));

    handles.RGB2(:,:,1) = 256 - R;
    handles.RGB2(:,:,2) = 256 - G;
    handles.RGB2(:,:,3) = 256 - B;

    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);
    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

function editSize_Callback(hObject, eventdata, handles)
% hObject    handle to editSize (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editSize as text
%        str2double(get(hObject,'String')) returns contents of editSize as a double

% --- Executes during object creation, after setting all properties.
function editSize_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editSize (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function editComment_Callback(hObject, eventdata, handles)
% hObject    handle to editComment (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editComment as text
%        str2double(get(hObject,'String')) returns contents of editComment as a double

% --- Executes during object creation, after setting all properties.
function editComment_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editComment (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function editPath_Callback(hObject, eventdata, handles)
% hObject    handle to editPath (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editPath as text
%        str2double(get(hObject,'String')) returns contents of editPath as a double

% --- Executes during object creation, after setting all properties.
function editPath_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editPath (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on slider movement.
function sliderBright_Callback(hObject, eventdata, handles)
% hObject    handle to sliderBright (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
set(handles.editBright,'String', sprintf('%10s:%4.0f%%', 'Brightness', 100*get(handles.sliderBright,'Value')));

% --- Executes during object creation, after setting all properties.
function sliderBright_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sliderBright (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function sliderContrast_Callback(hObject, eventdata, handles)
% hObject    handle to sliderContrast (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
set(handles.editContrast,'String', sprintf('%10s:%4.0f%%', 'Contrast', 100*get(handles.sliderContrast,'Value')));

% --- Executes during object creation, after setting all properties.
function sliderContrast_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sliderContrast (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function editBright_Callback(hObject, eventdata, handles)
% hObject    handle to editBright (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editBright as text
%        str2double(get(hObject,'String')) returns contents of editBright as a double

% --- Executes during object creation, after setting all properties.
function editBright_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editBright (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function editContrast_Callback(hObject, eventdata, handles)
% hObject    handle to editContrast (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editContrast as text
%        str2double(get(hObject,'String')) returns contents of editContrast as a double

% --- Executes during object creation, after setting all properties.
function editContrast_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editContrast (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbuttonContrastBrightness.
function pushbuttonContrastBrightness_Callback(hObject, eventdata, handles)
% hObject    handle to pushbuttonContrastBrightness (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if (handles.fileLoaded==1)
    handles.RGB2 = changeBrightness(handles.RGB, get(handles.sliderBright, 'Value'), get(handles.sliderContrast, 'Value'));
    axes(handles.axes2); imshow(handles.RGB2);
    handles.fileLoaded2 = 1;
    handles = updateHistograms(handles);
    guidata(hObject, handles);
else
    h = msgbox('No primary file has been loaded!','Error','error');
end

function handlesNew = updateHistograms(handles)
handlesNew = handles;
if (handles.fileLoaded == 1)
    set(handles.textHist1, 'Visible', 'on');
    axes(handlesNew.axesHist1);
    cla;
    ImageData1 = reshape(handlesNew.RGB(:,:,1), [size(handlesNew.RGB, 1) * size(handlesNew.RGB, 2) 1]);
    ImageData2 = reshape(handlesNew.RGB(:,:,2), [size(handlesNew.RGB, 1) * size(handlesNew.RGB, 2) 1]);
    ImageData3 = reshape(handlesNew.RGB(:,:,3), [size(handlesNew.RGB, 1) * size(handlesNew.RGB, 2) 1]);
    [H1, X1] = hist(ImageData1, 1:5:256);
    [H2, X2] = hist(ImageData2, 1:5:256);
    [H3, X3] = hist(ImageData3, 1:5:256);
    hold on;
    plot(X1, H1, 'r');
    plot(X2, H2, 'g');
    plot(X3, H3, 'b');
    axis([0 256 0 max([H1 H2 H3])]);
end
if (handlesNew.fileLoaded2 == 1)
    set(handles.textHist2, 'Visible', 'on');
    axes(handlesNew.axesHist2);
    cla;
    ImageData1 = reshape(handlesNew.RGB2(:,:,1), [size(handlesNew.RGB2, 1) * size(handlesNew.RGB2, 2) 1]);
    ImageData2 = reshape(handlesNew.RGB2(:,:,2), [size(handlesNew.RGB2, 1) * size(handlesNew.RGB2, 2) 1]);
    ImageData3 = reshape(handlesNew.RGB2(:,:,3), [size(handlesNew.RGB2, 1) * size(handlesNew.RGB2, 2) 1]);
    [H1, X1] = hist(ImageData1, 1:5:256);
    [H2, X2] = hist(ImageData2, 1:5:256);
    [H3, X3] = hist(ImageData3, 1:5:256);
    hold on;
    plot(X1, H1, 'r');
    plot(X2, H2, 'g');
    plot(X3, H3, 'b');
    axis([0 256 0 max([H1 H2 H3])]);
end

% --- Executes on button press in RotateButton.
function RotateButton_Callback(hObject, eventdata, handles)
% hObject    handle to RotateButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.RGB2 = imrotate(handles.RGB, round(get(handles.sliderRotate,'Value')));
handles.fileLoaded = 1;

axes(handles.axes2); cla; imshow(handles.RGB2);

handles = updateHistograms(handles);

guidata(hObject, handles);

% --- Executes on slider movement.
function sliderRotate_Callback(hObject, eventdata, handles)
% hObject    handle to sliderRotate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
set(handles.editRotate,'String', sprintf('%10s:%4.0f', 'Rotate', round(get(handles.sliderRotate,'Value'))));

% --- Executes during object creation, after setting all properties.
function sliderRotate_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sliderRotate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function editRotate_Callback(hObject, eventdata, handles)
% hObject    handle to editRotate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editRotate as text
%        str2double(get(hObject,'String')) returns contents of editRotate as a double

% --- Executes during object creation, after setting all properties.
function editRotate_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editRotate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


text:
function [dx] = spm_sde_dx(dfdx,dfdw,f,t)
% returns dx(t) = (expm(dfdx*t) - I)*inv(dfdx)*f + w for SDEs
% FORMAT [dx] = spm_sde_dx(dfdx,dfdw,f,t)
% dfdx   = df/dx - x: states
% dfdw   = df/dw - w: i.i.d. Weiner process
% f      = dx/dt
% t      = integration time: (default t = 1);
%
% dx     = x(t) - x(0)
%--------------------------------------------------------------------------
% Integration of stochastic differential equations using local linearization.
% This scheme accommodates nonlinearities in the state equation by using a
% functional of f(x) = dx/dt.  This uses the equality
%
%             expm([0    0]*t) = expm(dfdx*t) - I)*inv(dfdx)*f
%                  [f dfdx]
%
% When t -> Inf this reduces to
%
%              dx(t) = -inv(dfdx)*f
%
% for the SDE:  dx = dfdx*x*dt + sqrt(2)*dfdw*dw
%
% where w is a standard Wiener process. Unstable modes are removed using
% the systems eigenmodes.
%
% see also spm_dx
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_sde_dx.m 5932 2014-03-28 10:04:32Z karl $

% defaults
%--------------------------------------------------------------------------
if nargin < 3, t = 1; end

% compute stochastic terms {E = exp(dfdx*t), e = exp(dfdx*dt)}
%--------------------------------------------------------------------------
dfdx  = full(dfdx);
m     = length(dfdx);
N     = 256;
dt    = t/N;

% condition unstable eigenmodes
%--------------------------------------------------------------------------
[v,s] = eig(full(dfdx),'nobalance');
s     = diag(s);
u     = pinv(v);
s     = 1j*imag(s) + min(real(s),-4);
eJdt  = real(v*diag(exp(s*dt))*u);
dfdx  = real(v*diag(s)*u);

% flow operators
%--------------------------------------------------------------------------
eJt   = eJdt;
Q     = sparse(m,m);
R     = dfdw*dfdw'*2;
dQ    = eJt*R*eJt';
TOL   = trace(dQ)/64;
for i = 1:N

    % integrate and update exp(dfdx*t)
    %----------------------------------------------------------------------
    Q   = Q + dQ*dt;
    eJt = eJt*eJdt;
    dQ  = eJt*R*eJt';

    % convergence
    %----------------------------------------------------------------------
    if trace(dQ) < TOL, break, end

end

% scaled [Wiener] innovation
%--------------------------------------------------------------------------
w     = spm_sqrtm(Q)*randn(m,1);

% local linear solution plus stochastic term
%==========================================================================
dx    = spm_dx(dfdx,f,t) + w;

author: spm


text:
function [rgb,points3d,points3dMatrix]=read_3d_pts_general(depthInpaint,K,depthInpaintsize,imageName,crop)
    % K is [fx 0 cx; 0 fy cy; 0 0 1];
    % for uncrop image crop =[1,1];
    % imageName is the full path to image
    cx = K(1,3); cy = K(2,3);
    fx = K(1,1); fy = K(2,2);
    invalid = depthInpaint==0;
    if ~isempty(imageName)
        im = imread(imageName);
        rgb = im2double(im);
    else
        rgb =double(cat(3,zeros(depthInpaintsize(1),depthInpaintsize(2)),...
                    ones(depthInpaintsize(1),depthInpaintsize(2)),...
                    zeros(depthInpaintsize(1),depthInpaintsize(2))));
    end
    rgb = reshape(rgb, [], 3);
    %3D points
    [x,y] = meshgrid(1:depthInpaintsize(2), 1:depthInpaintsize(1));
    x3 = (x-cx).*depthInpaint*1/fx;
    y3 = (y-cy).*depthInpaint*1/fy;
    z3 = depthInpaint;
    points3dMatrix =cat(3,x3,z3,-y3);
    points3dMatrix(cat(3,invalid,invalid,invalid))=NaN;
    points3d = [x3(:) z3(:) -y3(:)];
    points3d(invalid(:),:) =NaN;
end

author: thusiyuan


text:
function tec_io_test02 ( )

%*****************************************************************************80
%
%% TEC_IO_TEST02 tests TEC_WRITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2006
%
%  Author:
%
%    John Burkardt
%
  tec_file_name = 'tiny.dat';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEC_IO_TEST02\n' );
  fprintf ( 1, '  TEC_WRITE can write finite element data to a TECPLOT ASCII file.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this example, we will write data to "%s".\n', tec_file_name );

  dim_num = 2;
  node_num = 5;
  element_num = 3;
  element_order = 3;
  node_data_num = 2;
  node_coord = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    2.0, 0.0; ...
    0.0, 1.0; ...
    1.0, 1.0 ]';
  element_node = [ ...
    1, 2, 4;
    5, 4, 2;
    2, 3, 5 ]';
  node_data = [ ...
    1.0, 0.0; ...
    0.8, 0.2; ...
    0.6, 0.4; ...
    0.9, 0.1; ...
    0.5, 0.5 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension         = %d\n', dim_num );
  fprintf ( 1, '  Number of nodes           = %d\n', node_num );
  fprintf ( 1, '  Number of elements        = %d\n', element_num );
  fprintf ( 1, '  Element order             = %d\n', element_order );
  fprintf ( 1, '  Number of node data items = %d\n', node_data_num );

  r8mat_transpose_print ( dim_num, node_num, node_coord, ...
    '  Coordinates of nodes:' );

  i4mat_transpose_print ( element_order, element_num, element_node, ...
    '  Nodes of elements:' );

  r8mat_transpose_print ( node_data_num, node_num, node_data, ...
    '  Node data for nodes:' );

  tec_write ( tec_file_name, dim_num, node_num, element_num, ...
    element_order, node_data_num, node_coord, element_node, node_data );

  return
end


text:
function vi = variablesPack(x,z)

d = size(x,2);
nx = size(x,1);
nz = size(z,1);

[x_indizes, z_indizes] = ocl.collocation.indizes(nx,nz,d);

vi = zeros(d*(nx+nz), 1);

vi(x_indizes) = x;
vi(z_indizes) = z;

end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function [fk,ier]=nufft2d3(nj,xj,yj,cj,iflag,eps,nk,sk,tk)
%NUFFT2D3: Nonuniform FFT in R^2 - Type 3.
%
%  [FK,IER] = NUFFT2D3(NJ,XJ,YJ,CJ,IFLAG,NK,SK,TK);
%
%                 1  nj
%     fk(k)    = -- SUM cj(j) exp(+/-i s(k) xj(j)) exp(+/-i t(k) yj(j))
%                nj j=1
%
%     If (iflag .ge.0) the + sign is used in the exponential.
%     If (iflag .lt.0) the - sign is used in the exponential.
%
%  Input parameters:
%
%     nj     number of sources   (integer)
%     xj,yj  location of sources (real *8)
%
%            on interval [-pi,pi].
%
%     cj     strengths of sources (complex *16)
%     iflag  determines sign of FFT (see above)
%     eps    precision request  (between 1.0e-15 and 1.0e-1)
%     nk     number of (noninteger) Fourier modes computed
%     sk,tk  k-values (locations) of desired Fourier modes
%
%  Output parameters:
%
%     fk     Fourier transform values (complex *16)
%     ier    error return code
%            ier = 0  => normal execution.
%            ier = 1  => precision eps requested is out of range.
%
%

fk=zeros(nk+3,1)+1i*zeros(nk+3,1);
ier=0;
% To avoid error with points aligned, add three points not aligned :
% [0 0], [0 1], [1 0]
sk = [sk;0;0;1];
tk = [tk;0;1;0];
xj = [xj;0;0;1];
yj = [yj;0;1;0];
nj = nj+3;
cj = [cj;0;0;0];
mex_id_ = 'nufft2d3f90(i int[x], i double[], i double[], i dcomplex[], i int[x], i double[x], i int[x], i double[], i double[], io dcomplex[], io int[x])';
[fk, ier] = nufft2d(mex_id_, nj, xj, yj, cj, iflag, eps, nk+3, sk, tk, fk, ier, 1, 1, 1, 1, 1);
fk = fk(1:end-3);


text:
function [W phi]=train_rbf(X,Y,Xc,k_i,basisfunction)
%trains a radial basis function
%X is a N_p by N_dim matrix of training data
%Y is a N_p by N_dim matrix of training data
%Xc is a N_r by N_dim matrix of rbf centres
%basisfunction may be 'gaussian' or 'polyharmonicspline'
%k_i is a prescaler for 'gaussian' rbf and function order for
%'polyharmonicspline'. See k_i(i)=0 for constant bias

%    Copyright Travis Wiens 2008
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%    If you would like to request that this software be licensed under a less
%    restrictive license (i.e. for commercial closed-source use) please
%    contact Travis at travis.mlfx@nutaksas.com

if nargin<4
    k_i=1;
end

if nargin<5
    basisfunction='gaussian';
end

N_r=size(Xc,1);%number of centres

W=zeros(N_r,1);%weight matrix
[z phi]=sim_rbf(Xc,X,W,k_i,basisfunction);%simulate rbf
A=pinv(phi'*phi)*phi';%do inverse
W=A*Y;%find weights

repo: mcbench-benchmarks