text:
function [vout, rout] = fitC13Data(v0, expdata, model, majorIterationLimit)
%
% USAGE:
%
%    [vout, rout] = fitC13Data(v0, expdata, model, majorIterationLimit)
%
% INPUTS:
%    v0:                     It will automatically be converted to alpha by solving `N*alpha = v`;
%                            if `v0` is a matrix then it is assumed to be a multiple start situation and
%                            `vout` will also have this size.
%    expdata:                either a data structure or a cell array of structures, in
%                            which case it is assumed that you wan to fit the sum of the scores
%    model:                  model structure
%
% OPTIONAL INPUT:
%    majorIterationLimit:    max number of iterations solver is allowed to take. Default = 1000
%
% OUTPUTS:
%    vout:                   reflects size of `v0`, result of NLPsolution
%    rout:                   cell, result of NLPsolution

if nargin < 4
    majorIterationLimit = 1000;
end
diffInterval = 1e-5;         %gradient step size.
method = 1; % 2 = in terms of v.  % 1 in terms of alpha
printLevel = 3; %3 prints every iteration.  1 does a summary.  0 = silent.

if method == 1
    if ~isfield(model, 'N')
       model.N = null(full(model.S));
       display('model.N should be defined');
       pause;
    end

    x0 = model.N\v0; % back substitute

    % safety check:
    if (max(abs(model.S*v0))> 1e-6)
        display('v0 not quite in null space');
        pause;
    end
    if(max(abs(model.N*x0 - v0)) > 1e-6)
        max(abs(model.N*x0 - v0))
        display('null basis is weird');
        pause;
    end

    % set up problem
    nalpha = size(model.N, 2);
    x_L = -1000*ones(nalpha,1);
    x_U = 1000*ones(nalpha,1);
    [A, b_L, b_U] = defineLinearConstraints(model, method);
elseif method == 2
    x0 = v0; % back substitute
    [A, x_L, x_U] = defineLinearConstraints(model, method);
    b_L = zeros(size(A,1),1);
    b_U = zeros(size(A,1),1);
else
    display('error'); pause;
end

numpoints = size(x0,2);
vout = zeros(size(v0));
rout = cell(numpoints, 1);

for k = 1:numpoints
    x_0 = x0(:,k);
    NLPproblem.objFunction = 'errorComputation2';
    NLPproblem.gradFunction = 'errorComputation2_grad';
    NLPproblem.lb = x_L;
    NLPproblem.ub = x_U;
    NLPproblem.name = 'c13fitting';
    NLPproblem.x0 = x_0;
    NLPproblem.A = A;
    NLPproblem.b_L = b_L;
    NLPproblem.b_U = b_U;
    NLPproblem.user.expdata = expdata;
    NLPproblem.user.model = model;
    NLPproblem.user.useparfor = true;
    NLPproblem.user.diff_interval = diffInterval;
    NLPproblem.osense = 1;

    NLPproblem.PriLevOpt = 1;
    cnan = ( method == 2);

    NLPsolution = solveCobraNLP(NLPproblem, 'checkNaN', cnan, 'printLevel', printLevel, 'iterationLimit', majorIterationLimit, 'logFile', 'minimize_SNOPT.txt');

    if method == 1
        vout(:,k) = model.N*NLPsolution.full;
    else
        vout(:,k) = NLPsolution.full;
    end
    rout{k} = NLPsolution;
end
return

author: opencobra

repo: cobratoolbox


text:
function RCM = DEM_demo_Bayesian_Model_Reduction
% This demonstration code illustrates the application of post hoc model
% optimisation or Bayesian model reduction (BMR) in identifying gene and
% gene-gene interaction effects in behavioural or physiological variables.
% The basic idea is to replace conventional heuristics based on the
% assumption that the contribution of any gene is sampled from a sparse
% distribution (such that a small number contribute and a large number do
% not) with an explicit search over a model space that includes all sparse
% and non-sparse models.  This exhaustive search rests upon recent
% advances in model optimisation based upon variational Bayesian model
% inversion.  In short, it is possible to estimate the posterior
% distribution of model parameters under a reduced model, given the
% posterior and prior distributions under a full model.
% In this context, a reduced model corresponds to a model in which some
% parameters are removed (by shrinking their prior variance to zero).
% This means that it is only necessary to invert the full model and then
% perform automatic BMR (over all possible combinations of parameters)
% using a greedy search based upon the free energy approximation to log
% model evidence.  With sufficient signal to noise, this scheme can
% recover the small number of effects, even in under determined or
% ill-posed problems (where the number of potential effects can vastly
% exceed the number of samples).
%
% The illustration below uses 128 subjects who have been measured three
% times (say in three brain regions) and we want to model these
% measurements in terms of first and second order genetic contributions
% given 8 (binary) genetic variables and all (unique) pair wise
% interactions.  This means that there are 36 unknown parameters
% (excluding a constant and, say, age confounds over subjects).  In the
% scheme below, each measurement is inverted separately under a simple
% (polynomial) model with uninformative priors on the parameters and
% (precision) hyper-parameters describing beliefs about signal to noise.
% A fixed effects Bayesian model averaging (BMA) scheme is used in
% combination with BMR to identify the best model out of all possible
% combinations of first and second order effects.  With the signal to
% noise and number of samples used in this simulation, the recovery is
% generally perfect.  This scheme also illustrates inference over a
% partition of model space (or families of models).

%__________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_Bayesian_Model_Reduction.m 6306 2015-01-18 20:50:38Z karl $

rng('default')

% Genomic data (G)
%--------------------------------------------------------------------------
Ns = 256;                                   % unmber of subjects
Ng = 8;                                    % unmber of genes
Nr = 3;                                    % unmber of (fMRI) regions

% Genotype
%--------------------------------------------------------------------------
G     = rand(Ns,Ng) > 1/2;

% design matix (U) with second-order (gene-gene interaction) effects
%--------------------------------------------------------------------------
GG    = [];
for i = 1:Ns
    u     = [];
    for j = 1:Ng
        for k = (j + 1):Ng
          u(end + 1) = G(i,j)*G(i,k);
        end
    end
    GG = [GG; u];
end
U  = [G GG];
Nb = size(U,2);
R  = ones(1,Nr);
A  = zeros(Nb,1);

% first-order effects (g1 abd g4) and one interaction (g1 x g2)
%--------------------------------------------------------------------------
A([1 4 (Ng + 3)],:) = [1/2 1 1];             % effect sizes

% Simulate data with a SNR of about 4:1
%--------------------------------------------------------------------------
y  = U*A*R + randn(Ns,Nr);                 % SD(noise) = 1

% Bayesian model inversion
%==========================================================================

% Model specification
%--------------------------------------------------------------------------
M.IS   = @(P,M,U) U*P.A;
M.pE.A = zeros(Nb,1);
M.pC.A = zeros(Nb,1) + 128;
M.hE   = 0;
M.hC   = 1/128;

% confounds (at the between subject level)
%--------------------------------------------------------------------------
age   = rand(Ns,1);
Y.X0  = [age ones(Ns,1)];

for i = 1:Nr

    % model inversion for this region
    %----------------------------------------------------------------------
    Y.y       = y(:,i);
    [Ep,Cp]   = spm_nlsi_GN(M,U,Y);

    % save model for subsequent BMS
    %----------------------------------------------------------------------
    DCM{i}.M  = M;
    DCM{i}.Ep = Ep;
    DCM{i}.Cp = Cp;

end

% Bayesian model reduction: see below for family fun
%==========================================================================
RCM  = spm_dcm_post_hoc(DCM,@fun);

% show results
%--------------------------------------------------------------------------
spm_figure('Getwin','Model posterior (over families)'); clf

Qf         = zeros(1,8);
Qf(fun(A)) = 1;                           % true family
Pf         = zeros(1,8);
Pf(1:length(RCM.Pf)) = RCM.Pf;            % psoteror family probailities

subplot(2,1,1)
bar([Pf; Qf]')
xlabel('familiy')
title('Model posterior (over families)','FontSize',16)
axis square
legend({'esimated','true'})

subplot(2,1,2)
bar([RCM.Ep.A A])
title('Recovered and true effects','FontSize',16)
xlabel('parameter')
axis square

% Partition of mdoel space into famllies
%==========================================================================
function k = fun(A)

f1  = any(A(1:4));
f2  = any(A(5:8));
f3  = any(A(9:12));

if ~f1 && ~f2 && ~f3, k = 1; end
if ~f1 && ~f2 &&  f3, k = 2; end
if ~f1 &&  f2 && ~f3, k = 3; end
if ~f1 &&  f2 &&  f3, k = 4; end
if  f1 && ~f2 && ~f3, k = 5; end
if  f1 && ~f2 &&  f3, k = 6; end
if  f1 &&  f2 && ~f3, k = 7; end
if  f1 &&  f2 &&  f3, k = 8; end

author: spm


text:
function [pd,Pdm,pmd] = som_probability_gmm(D, sM, K, P)

%SOM_PROBABILITY_GMM Probabilities based on a gaussian mixture model.
%
% [pd,Pdm,pmd] = som_probability_gmm(D, sM, K, P)
%
%   [K,P] = som_estimate_gmm(sM,D);
%   [pd,Pdm,pmd] = som_probability_gmm(D,sM,K,P);
%   som_show(sM,'color',pmd(:,1),'color',Pdm(:,1))
%
%  Input and output arguments:
%   D    (matrix) size dlen x dim, the data for which the
%        (struct) data struct,     probabilities are calculated
%   sM   (struct) map struct
%        (matrix) size munits x dim, the kernel centers
%   K    (matrix) size munits x dim, kernel width parameters
%                 computed by SOM_ESTIMATE_GMM
%   P    (matrix) size 1 x munits, a priori probabilities for each
%                 kernel computed by SOM_ESTIMATE_GMM
%
%   pd   (vector) size dlen x 1, probability of each data vector in
%                 terms of the whole gaussian mixture model
%   Pdm  (matrix) size munits x dlen, probability of each vector in
%                 terms of each kernel
%   pmd  (matrix) size munits x dlen, probability of each vector to
%                 have been generated by each kernel
%
% See also SOM_ESTIMATE_GMM.

% Contributed to SOM Toolbox vs2, February 2nd, 2000 by Esa Alhoniemi
% Copyright (c) by Esa Alhoniemi
% http://www.cis.hut.fi/projects/somtoolbox/

% ecco 180298 juuso 050100

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

% input arguments
if isstruct(sM), M = sM.codebook; else M = sM; end
[c dim] = size(M);

if isstruct(D), D = D.data; end
dlen = size(D,1);

% reserve space for output variables
pd = zeros(dlen,1);
if nargout>=2, Pdm = zeros(c,dlen); end
if nargout==3, pmd = zeros(c,dlen); end

% the parameters of each kernel
cCoeff = cell(c,1);
cCoinv = cell(c,1);
for m=1:c,
  co = diag(K(m,:));
  cCoinv{m} = inv(co);
  cCoeff{m} = 1 / ((2*pi)^(dim/2)*det(co)^.5);
end

% go through the vectors one by one
for i=1:dlen,

  x = D(i,:);

  % compute p(x|m)
  pxm = zeros(c,1);
  for m = 1:c,
    dx     = M(m,:) - x;
    pxm(m) = cCoeff{m} * exp(-.5 * dx * cCoinv{m} * dx');
    %pxm(m) = normal(dx, zeros(1,dim), diag(K(m,:)));
  end
  pxm(isnan(pxm(:))) = 0;

  % p(x|m)
  if nargin>=2, Pdm(:,i) = pxm; end

  % P(x) = P(x|M) = sum( P(m) * p(x|m) )
  pd(i) = P*pxm;

  % p(m|x) = p(x|m) * P(m) / P(x)
  if nargout==3, pmd(:,i) = (P' .* pxm) / pd(i); end

end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% subfunction normal
%
% computes probability of x when mean and covariance matrix
% of a distribution are known

function result = normal(x, mu, co)

[l dim] = size(x);
coinv   = inv(co);
coeff   = 1 / ((2*pi)^(dim/2)*det(co)^.5);
diff   = x - mu;
result = coeff * exp(-.5 * diff * coinv * diff');

author: ilarinieminen

repo: SOM-Toolbox


text:
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
c=@(x)(x.^2);
dc=@(x)(2*x);
Lx=1.0;
Nx=200;
dx=Lx/Nx;
m=createMesh1D(Nx, Lx);
x_face=m.facecenters.x;
x_cell=m.cellcenters.x;
BC=createBC(m);
BC.left.periodic=true;
BC.right.periodic=true;
u0=abs(sin(x_cell/Lx*10*pi));
u_old=createCellVariable(m, u0);
u_val=u_old;
dt=0.1;
c_face=createFaceVariable(m, 0.0);
c_face.xvalue=c(x_face);
dc_cell=createCellVariable(m, dc(x_cell));
Mconv=convectionUpwindTerm(c_face);
Ms=linearSourceTerm(dc_cell);
[Mbc, RHSbc]=boundaryCondition(BC);
for i=1:1000
  [Mt, RHSt]=transientTerm(u_old, dt, 1.0);
  M=Mt+Mconv-Ms+Mbc;
  RHS=RHSt+RHSbc;
  u_val=solvePDE(m, M, RHS);
  u_old=u_val;
  visualizeCells(u_val); drawnow;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [ m, d ] = easter_stewart ( y )

%*****************************************************************************80
%
%% EASTER_STEWART computes the month and day of Easter for a Gregorian year.
%
%  Example:
%
%    Y = 2001
%
%    A = 6
%    B = 20
%    C = 1
%    DD = 5
%    E = 0
%    G = 6
%    H = 18
%    MM = 0
%    J = 0
%    K = 1
%    L = 6
%    M = 4
%    D = 15
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Thomas O'Beirne,
%    Puzzles and Paradoxes,
%    Oxford University Press, 1965.
%
%    Ian Stewart,
%    Easter is a Quasicrystal,
%    Scientific American,
%    March 2001, pages 80-83.
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, integer M, D, the month and day of Easter.
%
  a = mod ( y, 19 );
  b = floor ( y / 100 );
  c = mod ( y, 100 );
  dd = floor ( b / 4 );
  e = mod ( b, 4 );
  g = floor ( ( 8 * b + 13 ) / 25 );
  h = mod ( 19 * a + b - dd - g + 15, 30 );
  mm = floor ( ( a + 11 * h ) / 319 );
  j = floor ( c / 4 );
  k = mod ( c, 4 );
  l = mod ( 2 * e + 2 * j - k - h + mm + 32, 7 );

  m = floor ( ( h - mm + l + 90 ) / 25 );
  d = mod ( h - mm + l + m + 19 , 32 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
filename='Cantileverbeam_Hexahedra_Linear_Structured_SYM';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';

nsteps = 10;
Vfrac_final = 0.15;
Perimeter_target=3.5;
optimality_final =1e-4;
constr_final =1e-4;

BCscale_factor = 0.3;
rotation_per_it = 20;
HJiter0 = 1;
e2 = 50;
N_holes = [8 4 4];
R_holes = 0.2;
phase_holes = [0 0 0];

Vfrac_initial = 0.7;
optimality_initial = 5e-2;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;

TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function repos = TestRepoList
    % TestRepoList. An example list of repositories for use with DepMatUpdate
    %
    %
    %
    %     Licence
    %     -------
    %     Part of DepMat. https://github.com/tomdoel/depmat
    %     Author: Tom Doel, 2015.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    repos = DepMatRepo.empty;
    repos(end + 1) = DepMatRepo('coremat', 'master', 'https://github.com/tomdoel/coremat.git', 'coremat_master');
    repos(end + 1) = DepMatRepo('dicomat', 'master', 'https://github.com/tomdoel/dicomat.git', 'dicomat_master');
    repos(end + 1) = DepMatRepo('matnat', 'master', 'https://github.com/tomdoel/matnat.git', 'matnat_master');
end

author: tomdoel

repo: pulmonarytoolkit


text:
clear;clc;
% caffe.set_mode_cpu();
caffe.set_mode_gpu();
caffe.set_device(0);
caffe.reset_all();

addpath('./evaluation_func/');
addpath('./evaluation_func/matlabPyrTools-master/');

% model = 'IDN_x2_deploy.prototxt';
% model = 'IDN_x4_deploy.prototxt';
model = 'IDN_x3_deploy.prototxt';

% weights = 'caffemodel/IDN_x2.caffemodel';
% weights = 'caffemodel/IDN_x4.caffemodel';
weights = 'caffemodel/IDN_x3.caffemodel';
%weights = 'caffemodel/IDN_x4_mscoco.caffemodel';
net=caffe.Net(model,weights,'test');
test_dataset='Set5'; % Set5 | Set14 | B100 | Urban100

testfolder=['test_data/' test_dataset '/'];
up_scale=3; % 2 | 3 | 4

savepath = 'results/';
folderResultCur = fullfile(savepath,[test_dataset,'_x',num2str(up_scale)]);
if ~exist(folderResultCur,'file')
    mkdir(folderResultCur);
end

if strcmp(test_dataset,'Set5') || strcmp(test_dataset,'Set14')
    filepaths=dir(fullfile(testfolder,'*.bmp'));
else
    filepaths=dir(fullfile(testfolder,'*.jpg'));
end

psnr_bic=zeros(length(filepaths),1);
psnr_idn=zeros(length(filepaths),1);

ssim_bic=zeros(length(filepaths),1);
ssim_idn=zeros(length(filepaths),1);

ifc_bic=zeros(length(filepaths),1);
ifc_idn=zeros(length(filepaths),1);

time_idn=zeros(length(filepaths),1);

for i=1:length(filepaths)
    %% read groud truth image
    [add,imname,type]=fileparts(filepaths(i).name);
    im=imread([testfolder imname type]);
    dimension=size(im,3);
    %% work on illuminance only
    if size(im,3)>1
        im_ycbcr=rgb2ycbcr(im);
        im=im_ycbcr(:,:,1);
        im_cb=im_ycbcr(:,:,2);
        im_cr=im_ycbcr(:,:,3);
        im_cb=modcrop(im_cb,up_scale);
        im_cr=modcrop(im_cr,up_scale);
        im_cb_=shave(imresize(imresize(im_cb,1/up_scale,'bicubic'),up_scale,'bicubic'),[up_scale,up_scale]);
        im_cr_=shave(imresize(imresize(im_cr,1/up_scale,'bicubic'),up_scale,'bicubic'),[up_scale,up_scale]);
    end
    im_gnd=modcrop(im,up_scale);
    im_gnd=single(im_gnd)/255;
    im_l=imresize(im_gnd,1/up_scale,'bicubic');

    %% bicubic interpolation
    im_b=imresize(im_l,up_scale,'bicubic');

    %% IDN feed-forward
    tic
    im_input=permute(im_l,[2,1,3]);
    net.blobs('data').reshape([size(im_input),1,1]);
    net.reshape();
    net.blobs('data').set_data(im_input);
    net.forward_prefilled();
    im_result=net.blobs('upsample').get_data();

    im_h=im_result'+mycrop(im_b,up_scale);
    time_idn(i)=toc;

    %% remove border
    im_h=myshave(uint8(im_h*255),up_scale);
    im_gnd=shave(uint8(im_gnd*255),[up_scale,up_scale]);
    im_b=shave(uint8(im_b*255),[up_scale,up_scale]);

    %% compute PSNR
    psnr_bic(i)=compute_psnr(im_gnd,im_b);
    psnr_idn(i)=compute_psnr(im_gnd,im_h);

    %% compute SSIM
    ssim_bic(i)=ssim_index(im_gnd,im_b);
    ssim_idn(i)=ssim_index(im_gnd,im_h);

    %% compute IFC

    ifc_idn(i) = ifcvec(double(im_gnd),double(im_h));
    ifc_bic(i)=ifcvec(double(im_gnd),double(im_b));

    %% save results
    if dimension>1
        imwrite(ycbcr2rgb(cat(3,im_h,im_cb_,im_cr_)),fullfile(folderResultCur,[imname,'_x',num2str(up_scale),'.png']));
    else
        imwrite(im_h,fullfile(folderResultCur,[imname,'_x',num2str(up_scale),'.png']));
    end
    save(fullfile(folderResultCur,['PSNR_',test_dataset,'_x',num2str(up_scale),'.mat']),'psnr_idn');
    save(fullfile(folderResultCur,['SSIM_',test_dataset,'_x',num2str(up_scale),'.mat']),'ssim_idn');
    save(fullfile(folderResultCur,['IFC_', test_dataset,'_x',num2str(up_scale),'.mat']),'ifc_idn');
end

fprintf('Mean PSNR for Bicubic: %f dB\n', mean(psnr_bic));
fprintf('Mean PSNR for IDN: %f dB\n', mean(psnr_idn));

fprintf('Mean SSIM for Bicubic: %f \n', mean(ssim_bic));
fprintf('Mean SSIM for IDN: %f \n', mean(ssim_idn));

fprintf('Mean IFC for Bicubic: %f \n', mean(ifc_bic));
fprintf('Mean IFC for IDN: %f \n', mean(ifc_idn));

fprintf('Mean Time for IDN: %f \n', mean(time_idn));

repo: IDN-Caffe


text:
function [signal, accept, typeeeg, rt, response, chan_names, pnts, nsweeps, rate, xmin, xmax]=loadeeg( FILENAME, chanlist, TrialList, typerange, acceptype, rtrange, responsetype)
% eeg_load_scan_eeg - Load Neuroscan .EEG format
%
% Usage: [signal, accept, typeeeg, rt, response, chan_names, pnts, nsweeps, rate, xmin, xmax]=loadeeg( FILENAME, chanlist, TrialList, typerange, acceptype, rtrange, responsetype)
%
%      FILENAME     input Neuroscan .avg file
%      signal	    output signal
%      variance     variance of the signal
%      chan_names   array that represent the name of the electrodes
%
%      i.e.
%	  [signal] = loadeeg( 'test.eeg' );     % load data into the array named 'signal'
%         plot( signal(1,:) );	  		% plot the signal for the first electrode of the first sweep
%
%      data are organised into an array of Number_of_electrode x (Number_of_points_per_trial*Number_of_sweeps)
%      for a file with 32 electrodes, 700 points per trial and 300 sweeps, the resulting array is
%      of 32 collumn and 700*300 rows (300 consecutive blocs of 700 points)
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:52 $

% Licence:  GNU GPL, no implied or express warranty
% History:  01/2001, arno_delorme@salk.edu
%	1062001		0.0		primitive version
%	1102001		1.0		fully working version
%	1112001		1.1		more parameters
%	1152001		1.2		fix bugs
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<1
	fprintf('Not enought arguments\n');
	help loadeeg
	return;
end;
if nargin<2 CHAN='all'; end;
if nargin<3 TrialList='all'; end;
if nargin<4 typerange='all'; end;
if nargin<5 acceptype='all'; end;
if nargin<6 rtrange  ='all'; end;
if nargin<7 responsetype='all'; end;

% open file for reading
% ---------------------
fid=fopen(FILENAME,'r','ieee-le');
if fid<0
	fprintf(2,['Error LOADEEG: File ' FILENAME ' not found\n']);
	return;
end;

% read general part of the erp header and set variables
% -----------------------------------------------------
erp = fread(fid, 362, 'uchar');	% skip the firsts 368 bytes
nsweeps = fread(fid, 1, 'ushort');	% number of sweeps
erp = fread(fid, 4, 'uchar'); 	% skip 4 bytes
pnts= fread(fid, 1, 'ushort');	% number of point per waveform
chan= fread(fid, 1, 'ushort');  % number of channels
erp = fread(fid, 4, 'uchar'); 	% skip 4 bytes
rate= fread(fid, 1, 'ushort');  % sample rate (Hz)
erp = fread(fid, 127, 'uchar');	% skip 125 bytes
xmin= fread(fid, 1, 'float32'); % in s
xmax= fread(fid, 1, 'float32'); % in s
erp = fread(fid, 387, 'uchar');	% skip 387 bytes

fprintf('number of channels         : %d\n', chan);
fprintf('number of points per trial : %d\n', pnts);
fprintf('sampling rate (Hz)         : %f\n', rate);
fprintf('xmin (s)                   : %f\n', xmin);
fprintf('xmax (s)                   : %f\n', xmax);

% read electrode configuration
% ----------------------------
fprintf('Electrode configuration\n');
for elec = 1:chan
   	channel_label_tmp = fread(fid, 10, 'uchar');
	chan_names(elec,:) = channel_label_tmp';
	for index = 2:9 if chan_names(elec,index) == 0 chan_names(elec,index)=' '; end; end;
	erp = fread(fid, 47-10, 'uchar');
	baseline(elec) = fread(fid, 1, 'ushort');
	erp = fread(fid, 10, 'uchar');
	sensitivity(elec) = fread(fid, 1, 'float32');
	erp = fread(fid, 8, 'uchar');
	calib(elec) = fread(fid, 1, 'float32');
	fprintf('%s: baseline: %d\tsensitivity: %f\tcalibration: %f\n', chan_names(elec,1:4), baseline(elec), sensitivity(elec), calib(elec));
	factor(elec) = calib(elec) * sensitivity(elec) / 204.8;
end;
%fprintf('Electrode configuration\n');
%for elec = 1:chan
%	erp = fread(fid, 47, 'uchar');
%	baseline(elec) = fread(fid, 1, 'ushort');
%	erp = fread(fid, 10, 'uchar');
%	sensitivity(elec) = fread(fid, 1, 'float32');
%	erp = fread(fid, 8, 'uchar');
%	calib(elec) = fread(fid, 1, 'float32');
%	fprintf('baseline: %d\tsensitivity: %f\tcalibration: %f\n', baseline(elec), sensitivity(elec), calib(elec));
%	factor(elec) = calib(elec) * sensitivity(elec) / 204.8;
%end;

xsize    = chan * pnts;
buf_size = chan * pnts ;			% size in shorts

% set tags for conditions
% -----------------------
if size(chanlist)  == size('all')	chanlist = [1:chan]; end;
if size(TrialList) == size('all')	trialtagI     = 1; else trialtagI     = 0; end;
if size(acceptype) == size('all')	acceptagI     = 1; else acceptagI     = 0; end;
if size(typerange) == size('all')	typetagI      = 1; else typetagI      = 0; end;
if size(responsetype) == size('all')	responsetagI  = 1; else responsetagI  = 0; end;
if size(rtrange)      == size('all')	rttagI        = 1; else rttagI        = 0; end;

count_selected = 1;
fprintf('Reserving array (can take some time)\n');
signal = zeros( chan, pnts*nsweeps);
fprintf('Array reserved, scanning file\n');

for sweep = 1:nsweeps

	% read sweeps header
	% ------------------
	s_accept   = fread(fid, 1, 'uchar');
	s_type     = fread(fid, 1, 'ushort');
	s_correct  = fread(fid, 1, 'ushort');
	s_rt       = fread(fid, 1, 'float32');
	s_response = fread(fid, 1, 'ushort');
	s_reserved = fread(fid, 1, 'ushort');

	unreaded_buf = 1;

	% store the sweep or reject the sweep
	% -----------------------------------
	if trialtagI trialtag = 1;        else trialtag = ismember(sweep, TrialList); end;
	if acceptagI acceptag = 1;        else acceptag =  ismember(s_accept, acceptype); end;
	if typetagI  typetag  = 1; 	  else typetag  =  ismember(s_type, typerange); end;
	if responsetagI responsetag  = 1; else responsetag  = ismember(s_response, responsetype); end;
	if rttagI       rttag  = 1; 	  else rttag  =  ismember(s_rt, rtrange); end;

	if typetag
		if trialtag
			if acceptag
				if responsetag
					if rttag

						buf = fread(fid, [chan pnts], 'short');
						unreaded_buf = 0;

						% copy information to array
						% -------------------------
						accept(count_selected)   = s_accept;
						typeeeg(count_selected)  = s_type;
						rt(count_selected)       = s_rt;
						response(count_selected) = s_response;

						% demultiplex the data buffer and convert to microvolts
						% -----------------------------------------------------
						for elec = 1:chan
							buf(elec, :) = (buf(elec, :)-baseline(elec)-0.0)*factor(elec);
						end;
						signal(:,[((count_selected-1)*pnts+1):count_selected*pnts]) = buf;
						count_selected = count_selected + 1;
						if not(mod(count_selected,10)) fprintf('%d sweeps selected out of %d\n', count_selected-1, sweep); end;
					end;
				end;
			end;
		end;
	end;

	if unreaded_buf fseek(fid, buf_size*2, 'cof'); end;
end;
nsweeps = count_selected-1;
fclose(fid);

% restrincting array
% ---------------------------------------
fprintf('rereservation of variables\n');
signal = signal(chanlist, 1:(count_selected-1)*pnts);
chan_names = chan_names(chanlist,:);

return;

% Frequency domain EEG File format
%
% The frequency domain epoched EEG format has
% the same header as that  described in the appendix
% of the SCAN manual (see the file sethead.h on the
% download page). For each sweep of data, there is a
% sweep header that is identical to that described
% on page Headers-5.
%
% At this point there is a difference.  After the sweep
% header, the frequency domain data for each sweep has
% the following format:
%
% for each channel (erp.nchannels):
%
% for each frequency bin (erp.pnts):
%
% real value of the FFT stored as a 4-byte float;
%
% for each frequency bin (erp.pnts):
%
% imaginary value of the FFT stored as a 4-byte float;
%
% The data have been scaled to microvolts prior to the FFT,
% and it is these results which are stored to the file.

author: PatternRecognition

repo: OpenBMI


text:
function p_homo = homogenize(p,y);
%HOMOGENIZE Homogenize polynomial
%
% f = homogenize(p,x)

deg   = degree(p);
deg_y = degree(y);
if rem(deg,deg_y)~=0
    error('The degree of the homogenizer is not an even fraction of deg(p).');
end

if 0
    error('The homogenizer must be homogenious.');
end

p_variables = getvariables(p);
p_homo = getbasematrix(p,0)*y^(deg/deg_y);
for i = 1:length(p_variables);
    monom = recover(p_variables(i));
    if degree(monom)<deg
        power = (deg-(degree(monom)))/deg_y;
        p_homo = p_homo + getbasematrix(p,p_variables(i))*monom*y^power;
    else
       p_homo = p_homo + getbasematrix(p,p_variables(i))*monom;
   end;
end
% Reset info about conic terms
p_homo.conicinfo = [0 0];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Gradient Boosting Regression 1D Example  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%               1) Load 1D Regression Datasets               %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (1a) (simple example) Generate Data from sine function
clear all; close all; clc;
nbSamples = 200;
epsilon   = 0.1;
x_limits  = [0, 100];

% Generate True function and data
X         = linspace(x_limits(1),x_limits(2),nbSamples)';
y_true    = sin(X*0.05);
y         = y_true + normrnd(0,epsilon,1,nbSamples)';

% Plot data
options             = [];
options.points_size = 10;
options.title       = 'noisy sinusoidal data';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1      = ml_plot_data([X(:),y(:)],options); hold on;

% Plot True function and Data
plot(X,y_true,'--k','LineWidth',2);
legend({'data','true function'})

%% (1b) (complex example) Generate Data from a sinc function
clear all; close all; clc;
% Set parameters for sinc function data
nbSamples = 200;
epsilon   = 0.1;
y_offset  = 0.5;
x_limits  = [-5, 5];

% Generate True function and data
X = linspace(x_limits(1),x_limits(2),nbSamples) ;
y_true = sinc(X) + y_offset ;
y = y_true + normrnd(0,epsilon,1,nbSamples);
X = X';
y = y';

% Plot data
options             = [];
options.points_size = 15;
options.title       = 'noisy sinc data';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1      = ml_plot_data([X(:),y(:)],options); hold on;

% Plot True function and Data
plot(X,y_true,'--k','LineWidth',2);
legend({'data','true function'})

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 2)  GRADIENT BOOSTING REGRESSION                     %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Gradient Boosting with decision stumps weak learners
% Options
options.nbWeakLearners = 200;

% Train Gradient Boosting Model
model = fitensemble(X,y,'LSBoost',options.nbWeakLearners,'Tree');

% Compute prediction

Y = predict(model,X);

% Plot data
options.title       = [options.title ' : Prediction'];

if exist('hl','var') && isvalid(hl), delete(hl);end
hl      = ml_plot_data([X(:),y(:)],options); hold on;

% Plot prediction from Gradient Boosting

plot(X,Y,'k','LineWidth',2);
legend({'data', 'predicted function'});

% Plot Loss function of Training Model

if exist('hl2','var') && isvalid(hl2), delete(hl2);end
figure
hl2 = plot(model.FitInfo);
xlabel('Number of weak learners','FontSize',14);
ylabel('Loss (MSE)','FontSize',14);
title('Cumulative Loss of the train data','FontSize',18);

author: epfl-lasa


text:
function a = ge( x, y )

evalin( 'caller', 'cvx_verify' );
b = cvx_pushcnstr( x, y, '>=' );
if nargout, a = b; end

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function movePTPEllipse_XY_1(iiwa,c,ratio,theta,velocity,accel,TefTool)
%% This funciton is used to
% Move end-effector of the robot on an ellipse the plane of the ellipse is
% parallel to the XY plane of the robot base.

%% Syntax:
% movePTPEllipse_XY(iiwa,c,ratio,theta,velocity,accel,TefTool)

%% Arreguments:
% iiwa: KST obejct.
% c: 2x1 vector,  is a vector defining the displacment of the the center
% point of the ellipse in the XY plane in relation to the start point of the ellipse,
% position units in (mm).
% ratio: the radious ratio (a/b) of the ellipse.
% theta: is the angle of the ellipe part , for drawing a complete
% ellipse this angle is equal to 2*pi.
% TefTool: 4X4 Transform matrix of the end-efector in the reference frame
% of the flange of the KUKA iiwa, position units in (mm).
% accel: is the acceleration (mm/sec2).
% velocity: is the motion velocity (mm/sec).

% Copy right: Mohammad SAFEEA
% 27th-June-2018

%% Convert inputs to a column vectors
c=colVec(c);
ratio=colVec(ratio);
theta=colVec(theta);
velocity=colVec(velocity);
accel=colVec(accel);

%% Check input variables
if(size(c,1)~=2)
fprintf('Error: the vector (c) shall be a 2x1 vector \n');
return;
end

if(size(ratio,1)~=1)
fprintf('Error: the ratio shall be a scalar  \n');
return;
end

if(size(theta,1)~=1)
fprintf('Error: the variable (theta) shall be a scalar  \n');
return;
end

if(size(TefTool,1)~=4)
fprintf('Error: the transform matrix (TefTool) shall be 4x4  \n');
return;
end

if(size(TefTool,2)~=4)
fprintf('Error: the transform matrix (TefTool) shall be 4x4  \n');
return;
end

if(size(velocity,1)~=1)
fprintf('Error: the velocity on the path shall be a scalar  \n');
return;
end

if(size(accel,1)~=1)
fprintf('Error: the acceleration on the path shall be a scalar  \n');
return;
end

dir=[1;0;0];
c=[c;
    0];

% perform the elliptical motion
movePTPEllipse_1(iiwa,c,dir,ratio,theta,velocity,accel,TefTool);
end

function y=colVec(x)
if size(x,2)==1
    y=x;
else
    y=x';
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function varargout = process_dipole_scanning( varargin )
% PROCESS_DIPOLE_SCANNING: Generates a brainstorm dipole file from the GLS and GLS-P inverse solutions.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Elizabeth Bock, John C. Mosher, Francois Tadel, 2013-2019

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Dipole scanning';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Sources';
    sProcess.Index       = 327;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/TutDipScan';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'results'};
    sProcess.OutputTypes = {'dipoles'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Definition of the options
     % === Time window
    sProcess.options.timewindow.Comment = 'Time window:';
    sProcess.options.timewindow.Type    = 'timewindow';
    sProcess.options.timewindow.Value   = [];
%     % === fit frequency
%     sProcess.options.downsample.Comment = 'Time between dipoles: ';
%     sProcess.options.downsample.Type    = 'value';
%     sProcess.options.downsample.Value   = {0.000, 'ms', []};
    % === Separator
    sProcess.options.sep2.Type = 'separator';
    sProcess.options.sep2.Comment = ' ';
    % === CLUSTERS
    sProcess.options.scouts.Comment = 'Limit scanning to selected scouts';
    sProcess.options.scouts.Type    = 'scout_confirm';
    sProcess.options.scouts.Value   = {};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInput) %#ok<DEFNU>
    OutputFiles = {};
    % === Get options
     % Time window to process
     %FitPeriod  = sProcess.options.downsample.Value{1};
    TimeWindow = sProcess.options.timewindow.Value{1};
    FitPeriod  = 0;
    AtlasList  = sProcess.options.scouts.Value;

    % === Get the sources
    % Read the source file
    sResultP = in_bst_results(sInput.FileName, 0);
    % Only accept dipole modeling in input
    switch sResultP.Function
        case {'gls_p','glsp','lcmvp'}
            % do nothing, valid
        otherwise
            bst_report('Error', sProcess, [], 'Dipole scanning is only available for the "Dipole modeling" and "LCMV" options (GLS-performance in the old interface).');
            return;
    end
    % Get the scouts structures
    if ~isempty(AtlasList)
        [sScouts, AtlasNames, sSurf] = process_extract_scout('GetScoutsInfo', sProcess, sInput, sResultP.SurfaceFile, AtlasList);
    else
        sScouts = [];
    end

    % === Get the results
    if isempty(sResultP.DataFile)
        DataMatP.Time = sResultP.Time;
        DataMatP.F = [];
        DataMatP.Comment = sResultP.Comment;
        % DataMatP.nAvg = sResultP.nAvg;
        DataMatP.Leff = sResultP.Leff;
    else
        DataMatP = in_bst_data(sResultP.DataFile);
        % Not supported for raw files
        if isstruct(DataMatP.F)
            bst_report('Error', sProcess, [], 'Dipole scanning is not supported for raw files. Import some data blocks first.');
            return;
        end
    end

    if isempty(sResultP.ImageGridAmp)
        % this can be massive, we should find a more iterative way
        % For now (March 2016), assume being only used for short data
        % sequences. (Mosher April 2016)
        sResultP.ImageGridAmp = sResultP.ImagingKernel * DataMatP.F(sResultP.GoodChannel,:);
    end

    % The function in_bst_results has already scaled the ImagingKernel to account for
    % data that has been averaged, so that the noise covariance scaling is correct.
    % So ImageGridAmp should also be correct.
    % But we need to calculate the possible Factor here for any other
    % manipulation of the data, as needed below.
    Factor = sqrt(DataMatP.Leff / sResultP.Leff); % is unity if both equal.
    if Factor ~= 1
        fprintf('BST Dipole Scanning> Need additional factor of %.2f (sqrt of %.1f) to account for average.\n',Factor, Factor^2);
    end

    % === Get the time
    if ~isempty(TimeWindow)
        SamplesBounds = bst_closest(TimeWindow, DataMatP.Time);
    else
        SamplesBounds = [1, size(sResultP.ImageGridAmp,2)];
    end
    timeVector = DataMatP.Time(SamplesBounds(1):SamplesBounds(2));

    % The "performance" image matrix
    Perf = sResultP.ImageGridAmp(:,SamplesBounds(1):SamplesBounds(2));
    % For constrained or mixed models this will flatten to norm
    sResultPFlat = process_source_flat('Compute',sResultP, 'rms');
    Pscan = sResultPFlat.ImageGridAmp(:,SamplesBounds(1):SamplesBounds(2));

    % === Find the index of 'best fit' at every time point
    % Get the selected scouts
    if ~isempty(sScouts)
        scoutVerts = [];
        for iScout = 1:length(sScouts)
            scoutVerts = [scoutVerts, sScouts(iScout).Vertices(:)'];
        end
        [mag,ind] = max(Pscan(scoutVerts,:));
        maxInd = scoutVerts(ind);
    else
        % don't use scouts, use all vertices
        [mag,maxInd] = max(Pscan,[],1);
    end

    % === Prepare a mask for downsampling the number of dipoles to save
    NumDipoles = size(Perf,2);
    if FitPeriod > 0
        tTotal = timeVector(end) - timeVector(1);
        nNewDipoles = tTotal/FitPeriod;
        dsFactor = round(NumDipoles/nNewDipoles);
    else
        dsFactor = 1;
    end
    temp = zeros(1,dsFactor);
    temp(1) = 1;
    dsMask = repmat(temp, 1, floor(NumDipoles/dsFactor));
    dsMask = logical([dsMask zeros(1,NumDipoles-length(dsMask))]);
    % downsample
    dsTimeVector = timeVector(dsMask);
    dsMaxInd = maxInd(dsMask);

    % So we have an index to each dipole location that has maximum
    % performance that we seek. Now find orientation and dipole amplitude

    % Get headmodel parameters
    sHeadModel = bst_get('HeadModelForStudy', sInput.iStudy);
    if isempty(sHeadModel)
        HeadModelFile = sResultP.HeadModelFile;
        if isempty(HeadModelFile)
            error('No headmodel available for this result file.');
        end
    else
        HeadModelFile = sHeadModel.FileName;
    end

    % HeadModelMat = in_bst_headmodel(HeadModelFile, 0, 'Gain', 'GridLoc', 'GridOrient');
    % Don't need large gain matrix
    HeadModelMat = in_bst_headmodel(HeadModelFile, 0, 'GridLoc', 'GridOrient');

    % === find the orientations
    switch (sResultP.nComponents)
        case 0  % MIXED HEAD MODEL
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
error('Not supported yet.');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %fullPerf = zeros(3, size(Perf,1), size(Perf,2));
            nTime = size(Perf,2);
            fullPerf = {};
            % Loop on all the regions
            for iScout = 1:length(sResultP.GridAtlas.Scouts)
                % Get the grid indices for this scouts
                iGrid = sResultP.GridAtlas.Scouts(iScout).GridRows;
                % If no vertices to read from this region: skip
                if isempty(iGrid)
                    continue;
                end
                % Get correpsonding row indices based on the type of region (constrained or unconstrained)
                switch (sResultP.GridAtlas.Scouts(iScout).Region(3))
                    case 'C'
                        fullPerf{end+1} = repmat(sResultP.GridOrient(iGrid,:)', 1, 1, nTime);
                    case {'U','L'}
                        % Convert from indices in GridLoc to indices in the source matrix
                        iSourceRows = bst_convert_indices(iGrid, sResultP.nComponents, sResultP.GridAtlas, 0);
                        fullPerf{end+1} = reshape(Perf(iSourceRows,:), 3, [], size(Perf,2));
                    otherwise
                        error(['Invalid region "' sResultP.GridAtlas.Scouts(iScout).Region '"']);
                end
            end
            % Concatenate the blocks
            fullPerf = cat(2, fullPerf{:});
            % Dipole orientation = amplitude in the three orientations
            orient = zeros(size(fullPerf,1),length(dsMaxInd)); % allocate
            for jj = 1:length(dsMaxInd)
                orient(:,jj) = fullPerf(:,dsMaxInd(jj),jj);
            end
        case 1
            % use the headmodel orientations
            kk = sResultP.nComponents; % notation convenience in below loop
            fullPerf = reshape(Perf,kk,[],size(Perf,2));
            orient = HeadModelMat.GridOrient(dsMaxInd,:)'; % simply the grid orient

            % now calculate the amplitude and embed in the orientation.
            % xhat = orientation * pinv(A * orientation) times data
            % A and data whitened and decomposed already, a rank one model
            for jj = 1:length(dsMaxInd)
                temp_ForwardModel = sResultP.SourceDecompSa(1,dsMaxInd(jj)); % scale factor from SVD
                temp_ForwardModel = temp_ForwardModel * Factor; % need to account for average vs epoched data.
                % Now pseudo invert this to get the scalar amplitude of source.
                % Note that fullPerf (the Perf value) already has the left singular
                % vectors.
                temp_estimate = fullPerf(1,dsMaxInd(jj),jj)/temp_ForwardModel;
                % replace unit orient with this orientation times the amplitude
                orient(:,jj) = temp_estimate * orient(:,jj);
            end
        case 2
            error('Not supported.');
        case 3
            % calculate the dipole orientations
            % first extract the calculation from the kernel, which is
            % performance:
            kk = sResultP.nComponents; % notation convenience in below loop
            fullPerf = reshape(Perf,kk,[],size(Perf,2));
            % Now calculate the optimal orientation from the performance
            orient = zeros(size(fullPerf,1),length(dsMaxInd)); % allocate
            for jj = 1:length(dsMaxInd)
                orient(:,jj) = fullPerf(:,dsMaxInd(jj),jj);
                % this is the raw orientation ("p") vector from the
                % performance measure. Since it is a single time slice (for
                % now April 2016), then no need to SVD. To get the optimal
                % orientation, we need to multiply it by the corresponding
                % head model decomposition components
                % Note, scale factor not important here since we normalize.
                orient(:,jj) = sResultP.SourceDecompVa(:,((1-kk):0)+(dsMaxInd(jj)*kk)) * diag(sResultP.SourceDecompSa(:,dsMaxInd(jj))) * orient(:,jj);
                orient(:,jj) = orient(:,jj) / norm(orient(:,jj)); % unit norm

                % with the orientation optimized, now calculate the
                % amplitude and embed in the orientation.
                % xhat = orientation * pinv(A * orientation) times data
                % A and data whitened and decomposed already, a rank three
                % model
                temp_ForwardModel = diag(sResultP.SourceDecompSa(:,dsMaxInd(jj))) * sResultP.SourceDecompVa(:,((1-kk):0)+(dsMaxInd(jj)*kk))'*orient(:,jj); % particular orientation
                temp_ForwardModel = temp_ForwardModel * Factor; % need to account for average vs epoched data.
                % Now pseudo invert this to get the scalar amplitude of source.
                % Note that fullPerf (the Perf value) already has the left singular
                % vectors.
                % pseudo invert temp_ForwardModel, a 3 x 1, by simply
                % transposing and dividing by the norm^2
                temp_estimate = (temp_ForwardModel'/ (temp_ForwardModel'*temp_ForwardModel)) * fullPerf(:,dsMaxInd(jj),jj);
                % replace unit orient with this orientation times the amplitude
                orient(:,jj) = temp_estimate * orient(:,jj);

            end
    end

%     % Dipole Amplitude, embed in orient. Use "fullPerf" generated above
%     %TBD: can be simplified from calculation given here, but this is direct
%     for jj = 1:length(dsMaxInd),
%         % form the source model, without the leading left singular vectors
%         kk = sResultP.nComponents; % notation convenience below
%         temp_ForwardModel = diag(sResultP.SourceDecompSa(:,dsMaxInd(jj))) * sResultP.SourceDecompVa(:,((1-kk):0)+(dsMaxInd(jj)*kk))'*orient(:,jj); % particular orientation
%         temp_ForwardModel = temp_ForwardModel * Factor; % need to account for average vs epoched data.
%         % Now pseudo invert this to get the scalar amplitude of source.
%         % Note that fullPerf (the Perf value) already has the left singular
%         % vectors.
%         temp_estimate = pinv(temp_ForwardModel) * fullPerf(:,dsMaxInd(jj),jj);
%         % replace unit orient with this orientation times the amplitude
%         orient(:,jj) = temp_estimate * orient(:,jj);
%     end
%
    % === Goodness of fit
    % square the performance at every source, therefore resulting in a
    % scalar squared performance value at every dipolar source
    %P2 = sum(reshape(abs(Perf).^2,sResultP.nComponents,[]),1);
    %P2 = reshape(P2,[],size(Perf,2));
    P2 = Pscan .^ 2;

    % get the squared norm of the whitened data, account for averaged data
    % factor as well.
    % Use the correct whitener, depending on lcmv vs gls
    switch sResultP.Options.InverseMethod
        case 'lcmv'
            wd2 = sum(abs(Factor * sResultP.DataWhitener * DataMatP.F(sResultP.GoodChannel,SamplesBounds(1):SamplesBounds(2))).^2,1);
        otherwise
            wd2 = sum(abs(Factor * sResultP.Whitener * DataMatP.F(sResultP.GoodChannel,SamplesBounds(1):SamplesBounds(2))).^2,1);
    end

    % the goodness of fit is now calculated by dividing the norm into each
    % performance value
    gof = P2 * diag(1./wd2);

    % === Chi-square
    % the chi square is the difference of the norm and the performance
    % resulting in the error for every source at every time point
    ChiSquare = repmat(wd2,size(P2,1),1) - P2;

    % The reduced chi-square is found by dividing by the degrees of freedom in
    % the error, which (for now) is simply a scalar, since we assume all
    % sources have the same degrees of freedom. Thus ROI modeling will require
    % that all ROIs have the same DOF.
    DOF = size(sResultP.ImagingKernel,2) - sResultP.nComponents;

    % downsample
    dsChiSquare = ChiSquare(:,dsMask);
    dsGOF = gof(:,dsMask);
    dsP = Pscan(:,dsMask);

    NumDipoles = length(dsMaxInd);

    %% === CREATE OUTPUT STRUCTURE ===
    bst_progress('start', 'Dipole File', 'Saving result...');
    % Get output study
    [sStudy, iStudy] = bst_get('Study', sInput.iStudy);
    % Comment: forced in the options
    if isfield(sProcess.options, 'Comment') && isfield(sProcess.options.Comment, 'Value') && ~isempty(sProcess.options.Comment.Value)
        Comment = sProcess.options.Comment.Value;
    % Comment: process default
    else
        Comment = [DataMatP.Comment ' | dipole-scan'];
    end
    % Get base filename
    [fPath, fBase, fExt] = bst_fileparts(sInput.FileName);
    % Create base structure
    DipolesMat = db_template('dipolemat');
    DipolesMat.Comment = Comment;
    DipolesMat.Time    = unique(dsTimeVector);

    % Fill structure
    for i = 1:NumDipoles
        DipolesMat.Dipole(i).Index          = 1;
        DipolesMat.Dipole(i).Time           = dsTimeVector(i);
        DipolesMat.Dipole(i).Origin         = [0 0 0];
        DipolesMat.Dipole(i).Loc            = HeadModelMat.GridLoc(dsMaxInd(i),:)';
        DipolesMat.Dipole(i).Amplitude      = orient(:,i);
        DipolesMat.Dipole(i).Errors         = 0;
        DipolesMat.Dipole(i).Noise          = [];
        DipolesMat.Dipole(i).SingleError    = [];
        DipolesMat.Dipole(i).ErrorMatrix    = [];
        DipolesMat.Dipole(i).ConfVol        = [];
        DipolesMat.Dipole(i).Probability    = [];
        DipolesMat.Dipole(i).NoiseEstimate  = [];
        DipolesMat.Dipole(i).Perform        = dsP(dsMaxInd(i),i);

        if ~isempty(dsGOF)
           DipolesMat.Dipole(i).Goodness = dsGOF(dsMaxInd(i),i);
        end
        if ~isempty(dsChiSquare)
           DipolesMat.Dipole(i).Khi2 = dsChiSquare(dsMaxInd(i),i);
        end
        if ~isempty(DOF)
           DipolesMat.Dipole(i).DOF = DOF;
        end
    end

    % Create the dipoles names list
    dipolesList = unique([DipolesMat.Dipole.Index]); %unique group names
    DipolesMat.DipoleNames = cell(1,length(dipolesList));
    k = 1; %index of names for groups with subsets
    nChanSet = 1;
    for i = 1:(length(dipolesList)/nChanSet)
        % If more than one channel subset, name the groups according to index and subset number
        if nChanSet > 1
            for j = 1:nChanSet
                DipolesMat.DipoleNames{k} = sprintf('Group #%d (%d)', dipolesList(i), j);
                DipolesMat.Subset(k) = j;
                k=k+1;
            end
        % If only one subsets, name the groups according to index
        else
            DipolesMat.DipoleNames{i} = sprintf('Group #%d', dipolesList(i));
            DipolesMat.Subset(i) = 1;
        end

    end
    DipolesMat.Subset = unique(DipolesMat.Subset);
    % Attach the new file to the input file
    DipolesMat.DataFile = sInput.FileName;
    % Add History field
    DipolesMat = bst_history('add', DipolesMat, 'scanning', ['Generated from: ' sInput.FileName]);

    % ===== SAVE NEW FILE =====
    % Create output filename
    ProtocolInfo = bst_get('ProtocolInfo');
    DipoleFile = bst_fullfile(ProtocolInfo.STUDIES, bst_fileparts(sStudy.FileName), 'dipoles_fit.mat');
    DipoleFile = file_unique(DipoleFile);
    % Save new file in Brainstorm format
    bst_save(DipoleFile, DipolesMat);

    % ===== UPDATE DATABASE =====
    % Create structure
    BstDipolesMat = db_template('Dipoles');
    BstDipolesMat.FileName = file_short(DipoleFile);
    BstDipolesMat.Comment  = Comment;
    BstDipolesMat.DataFile = sInput.FileName;
    % Add to study
    iDipole = length(sStudy.Dipoles) + 1;
    sStudy.Dipoles(iDipole) = BstDipolesMat;

    % Save study
    bst_set('Study', iStudy, sStudy);
    % Update tree
    panel_protocols('UpdateNode', 'Study', iStudy);
    % Save database
    db_save();
    % Return output file
    OutputFiles{1} = DipoleFile;
end

author: brainstorm-tools


text:
2     % problem
5     % outflow boundary
7     % grid parameter
3     % discretisation
0.02  % viscosity parameter
1     % Picard/Newton/hybrid linearization
9     % number of Picard iterations
1.d-5 % nonlinear tolerance

%% Data file for backward facing step


text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@sym uplus (@var{x})
%% @defopx Operator @@sym {+@var{x}} {}
%% Symbolic unitary minus.
%%
%% A no-op.  Example:
%% @example
%% @group
%% syms x
%% +x
%%   @result{} (sym) x
%% @end group
%% @end example
%% @seealso{@@sym/uminus}
%% @end defop

function x = uplus(x)

  % no-op

end

%!test
%! syms x
%! assert (isa (+x, 'sym'))
%! assert (isequal (+x, x))

%!test
%! A = sym([0 -1 inf]);
%! assert( isequal ( +A, A))

repo: octsympy


text:
%%
%
% This version is the working version from Pasadena 12/00
% it uses a single Mc cut, not a time variable cut.
% this is the most basic version of the pval code
%
%
%
%%

global m0 m events t dt obs_events calc_events maxt pl1
global bvm cut_cat all_minx fval tstep_tmp disc_events
global ncst length_log vts calc_time synth_events
global calc_sd  calc_ed pre_events daymc pl2  binv cua
global magco_fixed comp_flag

all_minx = [];
obs_events=[];
events = [];
synth_events = [];
cutobs_events = [];
nocutobs_events = [];
tstep_tmp = [];
vts = [];
daymc = [];
fval = [];
mmc = [];

report_this_filefun(mfilename('fullpath'));
m0 = maepi(1,6);
dt = .1;
%%
% get GUI input for time interval for calculation
% and for prediction.
%%

prompt   = {'Enter A,b,c,p calculation start time(days):',...
    'Enter calculation end time:',...
    'Enter forecast start time(days):',...
    'Enter forecast end time:'};
title    = 'Forecast based on A,b,c & p';
lines = 1;
def     = {'0','10','10','30'};
answer   = inputdlg(prompt,title,lines,def);
calc_sd = str2double(answer{1,1});
calc_ed = str2double(answer{2,1});
fore_sd = str2double(answer{3,1});
fore_ed = str2double(answer{4,1});

%%
% get catalog ffrom the plot cumulative number window
%%
tcut_cat = newt2;
hold_newt2 = newt2;
%%
% transform input days to decimal years
%%
calc_start = maepi(:,3) + calc_sd/365.0;
calc_end = maepi(:,3) + calc_ed/365.0;
fore_start = maepi(:,3) + fore_sd/365.0;
fore_end = maepi(:,3) + fore_ed/365.0;
endt = (calc_ed-calc_sd);

%%
% error check on GUI values
%%
maxt = (max(newt2.Date)-min(newt2.Date)+(0.05/365))*365;
if (calc_ed-calc_sd) > maxt
    disp('Less days in catalog that calculation period');
    return;
end

%%
% the following line is to take care of last eq's missed in
% bwith t calc
%%
m = min(newt2.Magnitude);

%%
% cut the events to the correct days
%
%%

ll = tcut_cat(:,3) <= calc_end;
cut_cat = tcut_cat(ll,:);
ll = cut_cat(:,3) < calc_start;
[pre_events,dump] = size(cut_cat(ll,3));
ll = cut_cat(:,3) > calc_start;
cut_cat=cut_cat(ll,:);

%%
% bvm from bwithtimc_nogui is used in the plot and c determination
% DO NOT USE THIS Mc
% calculate Mc in 100 event steps, returned in bvm
% same tool as for b with time
% defaults to best of 90 or 95% option
% ni is step size for Mc calc Set in bwithtimc_nogui too
%%
% use all events to compute mc
if ~exist('cutcat_flag', 'var')
    keep_newt2 = newt2;
    newt2 = cut_cat;
    bwithtimc_nogui;
    newt2 = keep_newt2;
end

%%
% calculate b value for entire sequence and get stand dev
% and use for constraints in search
%%
tmp_newt2=newt2;
newt2=cut_cat;
bdiff_omori;
bmin = bvml-.01*stanml;
bmax = bvml+.01*stanml;
bv = bvml;
bv0 = bvml;
av0 = avml;
newt2=tmp_newt2;
disp(['Initial b value: ',num2str(bvml),'+/-',num2str(stanml)]);

%%
% calculate a 6th order poly fit to the mc
% and find the mean of that for the last third
% of the data.  Used only for plotting
%%
[xb,yb] = size(bvm);
numpoly = min(xb-1,10);
csday = cumsum(bvm(:,5));
[pcof] = polyfit(csday,bvm(:,4),numpoly);
smoo_mc = polyval(pcof,csday);
mmc(1:length(csday)) = magco;

figure
pbvm = plot(csday,bvm(:,4),'-mo');
hold on
psm = plot(csday,smoo_mc,'-b');
pmag = plot(csday,mmc,'k-');
xlabel('days')
ylabel('Mc')
legend([pbvm,psm,pmag],'Mc in 100 event steps','Smoothed Mc','Mc cut','location','northeast');

tlen = calc_ed;

%%
% get GUI input for fixed mc
%%

prompt   = {'The calculated Mc is as follows, this value may be adjusted.'};
title    = 'Fixed Magnitude of Completeness';
lines = 1;
def     = {num2str(magco)};
answer   = inputdlg(prompt,title,lines,def);
magco_fixed = str2double(answer{1,1});

%%
% cut the catalog at the given Mc
%%
ll = cut_cat(:,6) >= magco_fixed - .05;
cut_cat = cut_cat(ll,:);

disp(['Number of events used: ',num2str(length(cut_cat))]);

%%
% see if cmin should be increased based on fixed Mc cutting
%%
cmin = .05;

for cib = 1:length(bvm)
    if bvm(cib,4) <= magco_fixed
        cmin = (bvm(cib,2)-maepi(1,3))*365
        break
    end
    if cib == length(bvm)
        disp('%%%ERROR%%% -- no bins complete to the Mc cut!!');
        return
    end
end

%create initial obs_events
binv = maepi(:,3)+0.05/365:dt/365:maepi(:,3)+calc_ed/365;
[nn,pdfbin] = histc(cut_cat(:,3),binv);
obs_events = cumsum(nn(1:length(nn)));

tbin = (binv-maepi(1,3))*365;
smc = interp1(csday,smoo_mc,tbin,'linear',smoo_mc(1));

%if compflag == 'of'
%    figure
%    pl1 = plot(obs_events,'-rx');
%    hold on
%    pl2 = plot(obs_events+50,'b');
%    legend([pl1,pl2],'Observed','Calculated Events',2);
%    drawnow
%end

vlb = [-3.5 bv0 cmin -0.30];
vub = [-1.0 bv0 cmin 1.99];
minx = [-2.00 bv0 cmin 2.0];

lobs = [];
Ap = [];
Aparmas = [];
Ap2=[];
lamda= [];
nlAp2=[];
obs = diff(obs_events(1:length(binv)));
dll = obs > 0;
nobs = obs(dll);
lobs = log10(nobs);

Mm=maepi(:,6)-magco_fixed;
cfixed = cmin;

stbin = tbin(2:length(tbin));
stbin = stbin(dll);
ll = length(stbin);
Aparams = zeros(ll,2);
Aparams(1:ll,1) = 1;
Aparams(1:ll,2) = Mm;

for i = 1:length(stbin)-1
    nlAp2(i) = ((stbin(i+1)-stbin(i)));
end
nlAp2 = [nlAp2 nlAp2(length(nlAp2'))]';
Ap2 = log10(nlAp2);

Ap = -log10(stbin+cfixed)';
Aparams = cat(2,Aparams,Ap);
id= eye(3);
id(1,1) = 1/10000;
id(3,3) = 1/10000;
id(2,2) = 1/10;
apri = [-1.67/10000 bv0/10 1.08/10000]';

lamda = (lobs-Ap2);
lamda = [lamda;apri];
Aparams = [Aparams; id];

valinv = Aparams\(lamda)

%options = optimset('Display','iter','Diagnostics','on','HessUpdate','dfp',...
%   'TolX',.000002,'TolFun',.0000001,'DiffMinChange',1e-3,...
%   'DiffMaxChange',.01)
%[x,fval] = fmincon('norm(func_omoris(x))',minx,[],[],[],[],vlb,vub,[],options);
%fval = [];

%rjA = log(10x(1))-bv0*(maepi(1,6)-magco_fixed);
x=[];
x(1) = valinv(1);
x(2) = bv0;
x(3) = cfixed;
x(4) = valinv(3);

if compflag == 'of'

    %%
    % plot results
    %%
    figure
    p1 = plot(newt2.Date,1:newt2.Count,'.k','MarkerSize',3);
    hold on
    l = cut_cat(:,3) > maepi(1,3) + x(3)/365;
    p2 = plot(cut_cat(l,3),1:length(cut_cat(l,3)),'.y','MarkerSize',6);
    p3 = plot(binv,obs_events,'r-.');
    %p4 = plot(binv,cumsum(calc_events),'b');

    obs = max(obs_events);

    forecast_omoris;
    p5 = plot(fore_time,(cumsum(fore_events)+obs+pre_events),'ms','Markersize',4);

    di = (max(obs_events) - sum(calc_events))

    legend([p1,p2,p3,p5],'Original','Cut in M','Observed','Forecast','Location', 'NorthWest');

    axes('Position',[.7,.1,.2,.2],'Box','on')
    axis off
    text(.1,.6,['A = ',num2str(x(1))]);
    text(.1,.45,['b = ',num2str(x(2))]);
    text(.1,.3,['c = ',num2str(x(3))]);
    text(.1,.15,['p = ',num2str(x(4))]);
end

newt2 = hold_newt2;

author: CelsoReyes


text:
im = iread('building2-1.png', 'grey');
idisp(im)
C = isift(im, 'nfeat', 200);

about C
C(1)
C(1).p

C(1).plot('y', 'scale', 16, 'clock')
idisp(im, 'dark')
C.plot('w', 'scale', 4, 'clock');

d = C(1).descriptor;
about d

author: petercorke

repo: machinevision-toolbox-matlab


text:
function varargout = auto_classify(varargin)
% AUTO_CLASSIFY M-file for auto_classify.fig
%      AUTO_CLASSIFY, by itself, creates a new AUTO_CLASSIFY or raises the existing
%      singleton*.
%
%      H = AUTO_CLASSIFY returns the handle to a new AUTO_CLASSIFY or the handle to
%      the existing singleton*.
%
%      AUTO_CLASSIFY('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in AUTO_CLASSIFY.M with the given input arguments.
%
%      AUTO_CLASSIFY('Property','Value',...) creates a new AUTO_CLASSIFY or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before auto_classify_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to auto_classify_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 auto_classify

% Last Modified by GUIDE v2.5 21-Jun-2006 00:34:13

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @auto_classify_OpeningFcn, ...
                   'gui_OutputFcn',  @auto_classify_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 auto_classify is made visible.
function auto_classify_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 auto_classify (see VARARGIN)

% input parameters

if isempty(varargin)
    handles.matrix2classify = [];
    handles.ncepstral = 5;
else
    handles.matrix2classify = varargin{1};
    handles.ncepstral = varargin{2};
end

handles.distancemeasure = 'sqEuclidean';

handles.cluster_method = 'kmeans';

set(handles.CepstralPopupMenu,'String',num2str([0:handles.ncepstral]'));
set(handles.CepstralPopupMenu,'Value',handles.ncepstral); %handles.ncepstral +1 );
% set(handles.CepstralPopupMenu,'Value',handles.ncepstral-1);

% Choose default command line output for auto_classify
% handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

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

uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = auto_classify_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;

close;

function clusters = clusterdata(X,k,dmeasure)
% A wrapper function for clustering

try
    clusters = kmeans(X,k,'distance',dmeasure,'EmptyAction','singleton','replicates',30);
catch
    clusters = k;
end

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

set(handles.AcceptButton,'Enable','off');
guidata(gcbo,handles);

cncepstral = get(handles.CepstralPopupMenu,'Value');

matrix2classify = handles.matrix2classify;
matrix2classify = madnormalize(matrix2classify,1); % normalize the columns
matrix2classify = matrix2classify(:,1:cncepstral+1); % include the cepstral coefficients that you want

if length(matrix2classify(1,:)) > 3
    matrix2classify = matrix2classify(:,[1,3:cncepstral+1]); % exclude the first cepestral coefficient
end

try
    segweight = str2num(get(handles.SegLenWeight,'String'));
catch
    segweight = 1;
    set(handles.SegLenWeight,'String','1');
end

% Weight the segment lengths by square root of the weight
matrix2classify(:,1) = matrix2classify(:,1)*sqrt(segweight);

% Other weightings can be added here but the cepestral coefficients have a
% natural exponential decline which acts a weighting

%Handle transformations to segment length

contents = get(handles.TransformPopupMenu,'String');
value = get(handles.TransformPopupMenu,'Value');

switch contents{value}
    case 'None'
        matrix2classify = matrix2classify;
    case 'Exclude'
        matrix2classify = matrix2classify(:,2:length(matrix2classify(1,:)));
    case 'Log'
        matrix2classify(:,1) = log(matrix2classify(:,1));
    case 'Exp'
        matrix2classify(:,1) = exp(matrix2classify(:,1));
end

if length(matrix2classify(1,:)) == 0 % matrix to classify has no information
    return
end

rangestate = get(handles.RangeSpecify,'Value');
diagnosticstate = get(handles.DiagnosticCheckbox,'Value');

if strcmp(handles.cluster_method,'hierarchical')
    handles.classification = cluster_hierarchical(matrix2classify);
    nclusters = length(unique(handles.classification)) - 1;
    set(handles.KClasses,'String',nclusters);

    if diagnosticstate
        pcaplot(matrix2classify,handles.classification);
    end

else if strcmp(handles.cluster_method, 'kmeans')

if not(rangestate)

    kclassstr = get(handles.KClasses,'String');

    if not(isempty(kclassstr))
        try
            nclasses = str2num(kclassstr);
        catch
            return;
        end

        classification =  clusterdata(matrix2classify,nclasses,handles.distancemeasure);% this can be changed to a different clustering algorithm

        % clusterdata shows how a clustering algorithm can be hooked in where the number of
        % classes need to be assigned.

        if length(classification) == 1
            msgbox('Error in clustering. Try changing the number of target clusters or rerunning.');
        else
            if diagnosticstate
                figure();
                [silh h] = silhouette(matrix2classify,classification);
                pcaplot(matrix2classify,classification);
            end
        handles.classification = classification;
        end
    else
        return;
    end
else
    try
        minclust = str2num(get(handles.MinClusters,'String'));
    catch
        return;
    end
    try
        maxclust = str2num(get(handles.MaxClusters,'String'));
    catch
        return;
    end

    resultsclust = cell((1+ maxclust)-minclust,1);

    hw = waitbar(0,'Clustering data set ');
    for i = 1:(1+maxclust)-minclust
        resultsclust{i} = clusterdata(matrix2classify,minclust + (i-1),handles.distancemeasure);
        if length(resultsclust{i}) == 1
            msgbox(['Clustering failed at ' num2str(resultsclust{i}) ' clusters. Try changing the range of clusters or rerunning.']);
            return;
        end
        waitbar(i/(1+maxclust-minclust));
    end
    close(hw);

    silhmean = zeros((1+ maxclust)-minclust,1);

    for i = 1:1+(maxclust)-minclust

        if diagnosticstate
            figure();
            [silh h] = silhouette(matrix2classify,resultsclust{i});
            pcaplot(matrix2classify,resultsclust{i});
        else
            silh = silhouette(matrix2classify,resultsclust{i});
        end

        silhmean(i) = mean(silh);
    end

    [x I] = max(silhmean);
    handles.classification = resultsclust{I(1)};

    set(handles.KClasses,'String',num2str(minclust + (I(1)-1)));
end
end
    end

set(handles.AcceptButton,'Enable','on');

guidata(gcbo,handles);

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

handles.output = 0;
guidata(hObject,handles);
uiresume(handles.figure1);

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

handles.output = handles.classification;
guidata(hObject,handles);
uiresume(handles.figure1);

function CepstralEdit_Callback(hObject, eventdata, handles)
% hObject    handle to CepstralEdit (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 CepstralEdit as text
%        str2double(get(hObject,'String')) returns contents of CepstralEdit as a double

% --- Executes during object creation, after setting all properties.
function CepstralEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to CepstralEdit (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function MinClusters_Callback(hObject, eventdata, handles)
% hObject    handle to MinClusters (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 MinClusters as text
%        str2double(get(hObject,'String')) returns contents of MinClusters as a double

% --- Executes during object creation, after setting all properties.
function MinClusters_CreateFcn(hObject, eventdata, handles)
% hObject    handle to MinClusters (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function MaxClusters_Callback(hObject, eventdata, handles)
% hObject    handle to MaxClusters (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 MaxClusters as text
%        str2double(get(hObject,'String')) returns contents of MaxClusters as a double

% --- Executes during object creation, after setting all properties.
function MaxClusters_CreateFcn(hObject, eventdata, handles)
% hObject    handle to MaxClusters (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in TransformPopupMenu.
function TransformPopupMenu_Callback(hObject, eventdata, handles)
% hObject    handle to TransformPopupMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns TransformPopupMenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from TransformPopupMenu

% value = get(hObject,'Value');
%
% if value == 2
%     handles.matrix2classify(:,1) = log(handles.matrix2classify(:,1));
% else
%     handles.matrix2classify(:,1) = exp(handles.matrix2classify(:,1));
% end

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

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

function KClasses_Callback(hObject, eventdata, handles)
% hObject    handle to KClasses (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 KClasses as text
%        str2double(get(hObject,'String')) returns contents of KClasses as a double

% --- Executes during object creation, after setting all properties.
function KClasses_CreateFcn(hObject, eventdata, handles)
% hObject    handle to KClasses (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in CepstralPopupMenu.
function CepstralPopupMenu_Callback(hObject, eventdata, handles)
% hObject    handle to CepstralPopupMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns CepstralPopupMenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from CepstralPopupMenu

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

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

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

% Hint: get(hObject,'Value') returns toggle state of RangeSpecify

state = get(hObject,'Value');

if state
    set(handles.KClasses,'Enable','off');
    set(handles.text3,'Visible','on');
    set(handles.MinClusters,'Visible','on');
    set(handles.text4,'Visible','on');
    set(handles.MaxClusters,'Visible','on');
else
    set(handles.KClasses,'Enable','on');
    set(handles.text3,'Visible','off');
    set(handles.MinClusters,'Visible','off');
    set(handles.text4,'Visible','off');
    set(handles.MaxClusters,'Visible','off');
end

% --- Executes on selection change in DistancePopupMenu.
function DistancePopupMenu_Callback(hObject, eventdata, handles)
% hObject    handle to DistancePopupMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns DistancePopupMenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from DistancePopupMenu

contents = get(hObject,'String');

dstring = contents{get(hObject,'Value')};

switch dstring
    case 'Squared Euclidean'
        handles.distancemeasure = 'sqEuclidean'
    case 'City Block (L1)'
        handles.distancemeasure = 'cityblock'
    case 'Cosine'
        handles.distancemeasure = 'cosine'
    case 'Correlation'
        handles.distancemeasure = 'correlation'
end

guidata(gcbo,handles);

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

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

handles.distancemeasure = 'sqEuclidean';

guidata(gcbo,handles);

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

% Hint: get(hObject,'Value') returns toggle state of DiagnosticCheckbox

function pcaplot(matrix2classify,classification)
% Plots the first two principal components uses different shapes selected
% randomly and different colors selected randomly to show classes.

nclusters = length(unique(classification));

if not(isempty(find(classification == 0)))
    nclusters = nclusters -1;
end

colors = rand(nclusters,3);
mark = mod(1:nclusters,13)+1;

[coef, score] = princomp(matrix2classify);

figure();

for i = 1:length(classification)
    h = line(score(i,1),score(i,2));
    markers = set(h,'Marker');
    if classification(i) == 0 % For segments which are not classified
        set(h,'Marker',markers{4});
        set(h,'Color',[0 0 0]);
    else
        set(h,'Marker',markers{mark(classification(i))});
        set(h,'Color',colors(classification(i),:));
    end
end

xlabel('PCA 1');
ylabel('PCA 2');

function classifymatrix = madnormalize(classifymatrix, cols2normalize)
  for i = cols2normalize
    classifymatrix(:,i) = (classifymatrix(:,i) - median(classifymatrix(:,i))) / mad(classifymatrix(:,i));
  end

% This code is modified from http://phys.columbia.edu/~aylin/clustering

function c=cluster_aylin(p1,p2);
%enter the p1,p2 found from hneighbors.m and get the struct array c.
%c(i).c will contain the indices of the objects that is in cluster i.
%c(1).c will be the cluster that has the maximum number of objects.
nn=size(p1,2);
p=std(p1);g=p(p2);p=p';
z=zeros(nn,1);zz=zeros(nn,1);
n=zeros(1,nn);
for i=1:nn
    n(i)=max(max(p1(:,p2(:,i))));
end
gg=[1:nn; max(g); n]';n=find(n<1);
[q1,q1]=sort(gg(n,2)');
g=gg(n(q1),1);
t=cputime;
j=1;m=1
	for i=1:length(n)
        ii=g(i);b=p2(:,ii);a=b;aa=a(find(gg(a,2)<1));
        if length([0; unique(z(aa))])<=2;a1=1;a2=0;
            while a1~=a2;a1=length(aa);
                a=unique(p2(:,a));a=a(find(gg(a,2)<1));
                a=unique(a(find(gg(a,2)<=mean(gg(aa,2))+std(gg(aa,2)))));
                a=a(find(ismember(a,aa)==0));
                if ~isempty(a);aa=[aa;a];
                    jj=aa(find(z(aa)));
                    if ~isempty(jj);;
                         u=unique(z(jj));
                         if length(u)==1;
                            zz(aa(find(~z(aa))))=m;z(aa)=u;m=m+1;
                         end
                         break;
                    end
                end;a2=length(aa);
            end;a=aa;
            jj=a(find(z(a)));
            if isempty(jj);
                z(a)=j;j=j+1;
                zz(a)=m;zz(b)=m-.1;zz(ii)=m-.2;m=m+1;
            end
        end
    end
    u=unique(z)
    u=u(find(u));v=length(u);vv=floor(1/v);if vv;vv=' is';else vv='s      are';end
    fprintf([int2str(v) ' cluster' vv ' found in ' num2str(cputime-t) 'sec\n\n'])
	q0=[];for i=1:v;qq=find(z==u(i));q0=[q0 length(qq)];end;
    [q1,q2]=sort(q0);%q2=q2(find(100*q1/nn>1));v=length(q2);
	c=[];for i=1:v;qq=find(z==u(q2(i)));[j,j]=sort(zz(qq));c(v-i+1).c=qq(j);end;
    if isempty(c);fprintf('\tNo cluster was found. \n \tscale the data (step size must be 1)\n');end

    function [p1,p2]=hneighbors(e);
	% this function finds the neighbors of each object in 'e' within a unit hypercube
	% and returns the sorted object distances to q.q1 and their identities to q.q2 , 'e' is a
	% matrix where i th row and j th column represents the j th component of the i th object.
	s='find(';i='abs(e(:,%d)-e(i,%d))<1&';
    % assuming that the data is given scaled and the characteristic step size is 1, variable s
	% keeps a script to find the objects that lie within a unit hypercube around the i th object.
    for j=1:size(e,2);
        s=[s sprintf(i,j,j)];
	end;s([end end+1])=');';
     % runs the script s for each of the objects and stores the sorted distances
     % from the i th object in q(i).q1 and their indentities in q(i).q2
    nn=size(e,1);m=ceil(nn^(1/4));p1=ones(m,nn);p2=kron(ones(m,1),1:nn);
    for i=1:nn;
        j=eval(s);
        [q,qq]=sort(sqrt(sum((e(j,:)-kron(ones(length(j),1),e(i,:)))'.^2)));q=q(find(q<1));mm=length(q);
        qq=j(qq(1:mm))';mn=min([m mm]);
        p1(1:mn,i)=q(1:mn)';
        p2(1:mn,i)=qq(1:mn)';
    end

function classification = cluster_hierarchical(matrix2classify)
    [p1,p2] = hneighbors(matrix2classify);
    c = cluster_aylin(p1,p2);

    nclusters = size(c,2);

    total_classified = 0

    nsegments = size(matrix2classify,1);

    classification = zeros(nsegments,1);

    for i = 1:nclusters
        for j = 1:length(c(i).c)
            classification(c(i).c(j)) = i;
        end
    end
;

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

% Hint: get(hObject,'Value') returns toggle state of HierarchicalRadioButton

selected = get(hObject,'Value');

if selected
    handles.cluster_method = 'hierarchical';
end

guidata(gcbo,handles);

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

% Hint: get(hObject,'Value') returns toggle state of KmeansRadioButton

selected = get(hObject,'Value');

if selected
    handles.cluster_method = 'kmeans';
end

guidata(gcbo,handles);

function SegLenWeight_Callback(hObject, eventdata, handles)
% hObject    handle to SegLenWeight (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 SegLenWeight as text
%        str2double(get(hObject,'String')) returns contents of SegLenWeight as a double

% --- Executes during object creation, after setting all properties.
function SegLenWeight_CreateFcn(hObject, eventdata, handles)
% hObject    handle to SegLenWeight (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [cm,cmd] = getColormaps

fnd = dir([mtex_path filesep 'plotting' filesep 'colormaps' filesep '*ColorMap.m']);

% remove extension to create command
cmd = cellfun(@(x) x(1:end-2),{fnd.name},'UniformOutput',false);

% remove ColorMap
cm = cellfun(@(x) x(1:end-8),cmd,'UniformOutput',false);

end

author: mtex-toolbox

repo: mtex

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


text:
function PIVlab_capture_chronos_save (cameraIP,nr_of_images,ImagePath,frame_nr_display)
hgui=getappdata(0,'hgui');
cameraURL = ['http://' cameraIP];
options = weboptions('MediaType','application/json','HeaderFields',{'Content-Type' 'application/json'});
resx=getappdata(hgui,'Chronos_resx');
resy=getappdata(hgui,'Chronos_resy');
bitdepth=getappdata(hgui,'Chronos_bits');
save_location=getappdata(hgui,'save_location');
save_type=getappdata(hgui,'save_type');

if matches(save_type,'TIFF')
	save_type_command='tiff';
elseif matches(save_type,'H264')
	save_type_command='h264';
elseif matches(save_type,'TIFF RAW')
	save_type_command='tiffraw';
end

%% save to SMB

if matches(save_location,'Download')
	% check if selected save method is existing:
	storage = webread ([cameraURL '/control/p/externalStorage']);
	if isfield(storage.externalStorage,'smb')
		set(frame_nr_display,'String',{'SMB File transfer setup OK!' ['Location: ' storage.externalStorage.smb.device]});drawnow;
		ok_go=1;
	else
		set(frame_nr_display,'String','Error: SMB File transfer not set up correctly');drawnow;
		ok_go=0;
	end
	if ok_go == 1
		if exist([storage.externalStorage.smb.device '\Chronos_PIVlab'],'dir') == 7
			[status,~] = rmdir([storage.externalStorage.smb.device '\Chronos_PIVlab'],'s');
			if status==0
				disp([storage.externalStorage.smb.device '\Chronos_PIVlab']')
				disp('Error: Temp directory could not be cleared.')
			end
		end

		response = webwrite([cameraURL '/control/startFilesave'],'filename','Chronos_PIVlab','format',save_type_command,'device','smb','start',1+2,'length',nr_of_images*2-1+2);

		pause(1)
		chronos_state = webread([cameraURL '/control/p/videoState']); %--> chronos_state.videoState 'filesave'
		%chronos_current_save_frame=webread([cameraURL '/control/p/playbackPosition']); % chronos_current_save_frame.playbackPosition

		exit_cntr=0;
		saving_running=1;
		while strcmp(chronos_state.videoState, 'recording') || strcmp(chronos_state.videoState, 'idle' ) || strcmp(chronos_state.videoState, 'live' )|| strcmp(chronos_state.videoState, 'paused')
			chronos_state = webread([cameraURL '/control/p/videoState']);
			exit_cntr=exit_cntr+1;
			pause(0.5)
			if exit_cntr >= 20
				disp('ERROR: timeout waiting for file transfer begin.')
				saving_running=0;
				response=webread([cameraURL '/control/stopFilesave']);
				break
			end
		end
		set(frame_nr_display,'String','Initiating transfer...');drawnow;
		msg_loop=0;
		if saving_running
			start_save=tic;
			while strcmp(chronos_state.videoState, 'filesave')
				if getappdata(hgui,'cancel_capture') ==1
					response=webread([cameraURL '/control/stopFilesave']);
					disp('cancelling save')
					break
				end
				pause(3)
				chronos_state = webread([cameraURL '/control/p/videoState']); %--> chronos_state.videoState 'filesave'
				direc=dir([storage.externalStorage.smb.device '\Chronos_PIVlab\*.tiff']);
				progress=min ([(length(direc) / (nr_of_images*2)*100) 99]);
				zeit=toc(start_save);
				done=progress/100;
				tocome=1-done;
				zeit=zeit/done*tocome;
				hrs=zeit/60^2;
				mins=(hrs-floor(hrs))*60;
				secs=(mins-floor(mins))*60;
				hrs=floor(hrs);
				mins=floor(mins);
				secs=floor(secs);
				first_line=['Download progress: ' num2str(round(progress*10)/10)  ' %'];
				if msg_loop>=2
					second_line=['Time left: ' sprintf('%2.2d', hrs) 'h ' sprintf('%2.2d', mins) 'm ' sprintf('%2.2d', secs) 's'];
				else
					second_line='...';
				end
				set(frame_nr_display,'String',{first_line second_line});drawnow;
				msg_loop=msg_loop+1;
			end
		end

		[status,~,~] = movefile([storage.externalStorage.smb.device '\Chronos_PIVlab\*.*'],ImagePath,"f");
		if status==0
			disp(['Temp dir: ' storage.externalStorage.smb.device '\Chronos_PIVlab']')
			disp(['Target dir: ' ImagePath])
			disp('Error: Could not move files from temp dir to target dir.')
		end
		[status,~] = rmdir([storage.externalStorage.smb.device '\Chronos_PIVlab'],'s');
		if status==0
			disp([storage.externalStorage.smb.device '\Chronos_PIVlab']')
			disp('Error: Temp directory could not be cleared.')
		end

		elapsed_time=toc(start_save);
		if elapsed_time < 121
			disp(['Elapsed time: ' num2str(round(elapsed_time)) ' seconds.'])
		else
			disp(['Elapsed time: ' num2str(round(elapsed_time/6)/10) ' minutes.'])
		end

		response=webread([cameraURL '/control/startLivedisplay']);
	end
	%% SAVE TO SSD or SD-Card

elseif matches(save_location,'SSD') || matches(save_location,'SD Card')

	if matches(save_location,'SSD')
		storage = webread ([cameraURL '/control/p/externalStorage']);
		if isfield(storage.externalStorage,'sda1')
			set(frame_nr_display,'String','SSD File transfer setup OK!');drawnow;
			ok_go=1;
			medium_string='SSD write progress: ';
		else
			set(frame_nr_display,'String','Error: SSD File transfer not set up correctly');drawnow;
			ok_go=0;
		end
	end

	if matches(save_location,'SD Card')
		storage = webread ([cameraURL '/control/p/externalStorage']);
		if isfield(storage.externalStorage,'mmcblk1p1')
			set(frame_nr_display,'String','SD-card file transfer setup OK!');drawnow;
			ok_go=1;
			medium_string='SD card write progress: ';
		else
			set(frame_nr_display,'String','Error: SD-card File transfer not set up correctly');drawnow;
			ok_go=0;
		end
	end

	if ok_go == 1
		if ImagePath(end) == filesep %remove trailing slashes if existent
			ImagePath(end) = [];
		end
		dindex=strfind(ImagePath,filesep);
		project_name=[ImagePath((dindex(end)+1:end)) '_' datestr(datetime(now,'ConvertFrom','datenum'),'YYYY_dd_mm_HH_MM_SS')];

		if matches(save_location,'SSD')
			response = webwrite([cameraURL '/control/startFilesave'],'filename',project_name,'format',save_type_command,'device','sda1','start',1+2,'length',nr_of_images*2-1+2);
		end
		if matches(save_location,'SD Card')
			response = webwrite([cameraURL '/control/startFilesave'],'filename',project_name,'format',save_type_command,'device','mmcblk1p1','start',1+2,'length',nr_of_images*2-1+2);
		end

		pause(1)
		chronos_state = webread([cameraURL '/control/p/videoState']); %--> chronos_state.videoState 'filesave'

		exit_cntr=0;
		saving_running=1;
		while strcmp(chronos_state.videoState, 'recording') || strcmp(chronos_state.videoState, 'idle' ) || strcmp(chronos_state.videoState, 'live' )|| strcmp(chronos_state.videoState, 'paused')
			chronos_state = webread([cameraURL '/control/p/videoState']);
			exit_cntr=exit_cntr+1;
			pause(0.5)
			if exit_cntr >= 20
				disp('ERROR: timeout waiting for file save begin.')
				saving_running=0;
				response=webread([cameraURL '/control/stopFilesave']);
				break
			end
		end
		set(frame_nr_display,'String','Initiating transfer...');drawnow;
		response = webwrite([cameraURL '/control/p'],'playbackPosition',0); %set to first frame, maybe this fixes jumping progress report...? No... :-(
		urlread([cameraURL '/control/p/totalFrames']); %fix jumping progress...? No...

		msg_loop=0;
		if saving_running
			start_save=tic;
			while strcmp(chronos_state.videoState, 'filesave')
				if getappdata(hgui,'cancel_capture') ==1
					response=webread([cameraURL '/control/stopFilesave']);
					disp('cancelling save')
					break
				end
				chronos_current_save_frame=webread([cameraURL '/control/p/playbackPosition']); % chronos_current_save_frame.playbackPosition
				chronos_state = webread([cameraURL '/control/p/videoState']); %--> chronos_state.videoState 'filesave'
				progress=min ([(chronos_current_save_frame.playbackPosition / (nr_of_images*2)*100) 99]);

				zeit=toc(start_save);
				done=progress/100;
				tocome=1-done;
				zeit=zeit/done*tocome;
				hrs=zeit/60^2;
				mins=(hrs-floor(hrs))*60;
				secs=(mins-floor(mins))*60;
				hrs=floor(hrs);
				mins=floor(mins);
				secs=floor(secs);
				first_line=[medium_string num2str(round(progress*10)/10)  ' %'];
				if msg_loop>=2
					second_line=['Time left: ' sprintf('%2.2d', hrs) 'h ' sprintf('%2.2d', mins) 'm ' sprintf('%2.2d', secs) 's'];
				else
					second_line='...';
				end
				set(frame_nr_display,'String',{first_line second_line});drawnow;
				msg_loop=msg_loop+1;
				pause(2);
			end
		end

	end
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function [nll,g,H,T] = LogisticLoss(w,X,y)
% w(feature,1)
% X(instance,feature)
% y(instance,1)

[n,p] = size(X);

Xw = X*w;
yXw = y.*Xw;

nll = sum(mylogsumexp([zeros(n,1) -yXw]));

if nargout > 1
    if nargout > 2
        sig = 1./(1+exp(-yXw));
        g = -X.'*(y.*(1-sig));
    else
        %g = -X.'*(y./(1+exp(yXw)));
        g = -(X.'*(y./(1+exp(yXw))));
    end
end

if nargout > 2
    H = X.'*diag(sparse(sig.*(1-sig)))*X;
end

if nargout > 3
    T = zeros(p,p,p);
    for j1 = 1:p
        for j2 = 1:p
            for j3 = 1:p
                T(j1,j2,j3) = sum(y(:).^3.*X(:,j1).*X(:,j2).*X(:,j3).*sig.*(1-sig).*(1-2*sig));
            end
        end
    end
end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function geometry_test0386 ( )

%*****************************************************************************80
%
%% TEST0386 tests LINES_EXP_EQUAL_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  test_num = 6;

  p1_test = [ ...
    0.0, 0.0; ...
    0.0, 0.0; ...
    0.0, 0.0; ...
    0.0, 0.0; ...
    0.0, 0.0; ...
    0.0, 0.0 ]';
  p2_test = [ ...
    1.0, 2.0; ...
    1.0, 2.0; ...
    1.0, 2.0; ...
    1.0, 2.0; ...
    1.0, 2.0; ...
    1.0, 2.0 ]';
  q1_test = [ ...
    0.0,  0.0; ...
    1.0,  2.0; ...
    0.0,  0.0; ...
    7.0, 14.0; ...
    1.0,  2.0; ...
    0.0, 10.0 ]';
  q2_test = [ ...
    1.0,  2.0; ...
    0.0,  0.0; ...
    2.0,  4.0; ...
    5.5, 11.0; ...
    3.0,  5.0; ...
    1.0, 12.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0386\n' );
  fprintf ( 1, '  LINES_EXP_EQUAL_2D tries to determine if two\n' );
  fprintf ( 1, '    explicit lines in 2D are equal.\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    p1(1:dim_num) = p1_test(1:dim_num,test)';
    p2(1:dim_num) = p2_test(1:dim_num,test)';
    q1(1:dim_num) = q1_test(1:dim_num,test)';
    q2(1:dim_num) = q2_test(1:dim_num,test)';

    fprintf ( 1, '\n' );
    fprintf ( 1, '  P1  %8f  %8f\n', p1(1:dim_num) );
    fprintf ( 1, '  P2  %8f  %8f\n', p2(1:dim_num) );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Q1  %8f  %8f\n', q1(1:dim_num) );
    fprintf ( 1, '  Q2  %8f  %8f\n', q2(1:dim_num) );

    equal = lines_exp_equal_2d ( p1, p2, q1, q2 );

    if ( equal )
      fprintf ( 1, '  The lines are equal.\n' );
    else
      fprintf ( 1, '  The lines are distinct.\n' );
    end

  end

  return
end


text:
function varargout = colsubset(varargin)
% VL_COLSUBSET Select a given number of columns
%   Y = VL_COLSUBSET(X, N) returns a random subset Y of N columns of
%   X. The selection is order-preserving and without replacement. If N
%   is larger or equal to the number of columns of X (e.g. N = Inf),
%   then the function returns all the columns (i.e., Y = X).
%
%   If 0 < N < 1, then the function returns a fraction N of the
%   columns (rounded to the closest integer).
%
%   [Y, SEL] = VL_COLSUBSET(...) returns the indexes SEL of the
%   selected columns.
%
%   The function accepts the following options:
%
%   Beginning::
%     Returns the fist N columns.
%
%   Ending::
%     Returns the last N columns.
%
%   Random:: [default]
%     Returns N columns selected at random (using RANDPERM()).
%
%   Uniform::
%     Returns N uniformly spaced columns.
%
%   Largest::
%     Returns the N largest columns (using SORTROWS()).
%
%   Smallest::
%     Returns the N smallest columns (using SORTROWS()).
%
%  See also: VL_HELP().
[varargout{1:nargout}] = vl_colsubset(varargin{:});

author: yihui-he

repo: panorama

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


text:
function outsig=rampdown(L,wintype)
%RAMPDOWN  Falling ramp function
%   Usage: outsig=rampdown(siglen);
%
%   `rampdown(siglen)` will return a falling ramp function of length
%   *siglen*. The ramp is a sinusoid starting from one and ending at
%   zero. The ramp is centered such that the first element is always one and
%   the last element is not quite zero, such that the ramp fits with
%   following zeros.
%
%   `rampdown(L,wintype)` will use another window for ramping. This may be
%   any of the window types from |firwin|. Please see the help on |firwin|
%   for more information. The default is to use a piece of the Hann window.
%
%   See also: rampup, rampsignal, firwin

complainif_argnonotinrange(nargin,1,2,mfilename);

if nargin==1
  wintype='hann';
end;

win=firwin(wintype,2*L,'inf');
outsig=win(1:L);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat