text:
function [ fea, out ] = ex_natural_convection( varargin )
%EX_NATURAL_CONVECTION 2D Example for natural convection of air in a square cavity.
%
%   [ FEA, OUT ] = EX_NATURAL_CONVECTION( VARARGIN ) Sets up and solves a natural convection
%   benchmark problem. Reference solutions are for example reported in G. Davis
%   "Natural convection of air in a square cavity a bench mark numerical solution",
%   IJNMF vol. 3, 249-254 (1983).
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       Ra          scalar {1e3}           Rayleigh number
%       Pr          scalar {0.71}          Prandtl number
%       l           scalar {1}             Side length of cavity
%       igrid       scalar 0/{1}           Cell type (0=quadrilaterals, 1=triangles)
%       hmax        scalar {0.05}          Max grid cell size
%       sf_u        string {sflag2}        Shape function for velocity
%       sf_p        string {sflag1}        Shape function for pressure
%       sf_T        string {sflag2}        Shape function for temperature
%       iphys       scalar 0/{1}           Use physics mode to define problem (=1)
%       iplot       scalar 0/{1}           Plot solution and error (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'Ra',       1e3; ...
  'Pr',       0.71; ...
  'l',        1; ...
  'igrid',    0; ...
  'hmax',     0.05; ...
  'sf_u',     'sflag2'; ...
  'sf_p',     'sflag1'; ...
  'sf_T',     'sflag2'; ...
  'iphys',    1; ...
  'iplot',    1; ...
  'tol',      0.1; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Reference values.
Ra_ref      = [ 1e3   1e4    1e5     1e6     ];
u_max_ref   = [ 3.649 16.178 34.73   64.63   ];
y_max_ref   = [ 0.813  0.823  0.855   0.850  ];
v_max_ref   = [ 3.697 19.617 68.59  219.36   ];
x_max_ref   = [ 0.178  0.119 0.066    0.0379 ];
Nu_mean_ref = [ 1.118  2.243 4.519    8.800  ];

% Model parameters.
Ra    = opt.Ra;      % Rayleigh number.
Pr    = opt.Pr;      % Prandtl number.
l     = opt.l;       % Length of rectangular domain.
sf_u  = opt.sf_u;    % FEM shape function type for velocity.
sf_p  = opt.sf_p;    % FEM shape function type for pressure.
sf_T  = opt.sf_T;    % FEM shape function type for temperature.

% Geometry definition.
fea.sdim         = { 'x' 'y' };
fea.geom.objects = { gobj_rectangle( 0, l, 0, l ) };

% Grid generation.
if( opt.igrid<=0 )
  fea.grid = rectgrid( round(l/opt.hmax), round(l/opt.hmax), [0 l;0 l] );
  if( opt.igrid<0 )
    fea.grid = quad2tri(fea.grid);
  end
else
  fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid );
end

% Boundary conditions.
n_bdr = max(fea.grid.b(3,:));    % Increment number of boundaries.
dtol  = opt.hmax;
ib_l  = findbdr( fea, ['x<',num2str(dtol)] );     % Right boundary number.
ib_r  = findbdr( fea, ['x>',num2str(l-dtol)] );   % Left boundary number.
ib_b  = findbdr( fea, ['y<',num2str(dtol)] );     % Bottom boundary number.
ib_t  = findbdr( fea, ['y>',num2str(l-dtol)] );   % Top boundary number.

% Add pressure point constraint on point closest to origin.
[~,ix] = min( fea.grid.p(1,:).^2 + fea.grid.p(2,:).^2 );
fea.pnt.index = ix;
fea.pnt.type  = 'constr';
fea.pnt.dvar  = 'p';
fea.pnt.expr  = 0';

% Problem definition.
if ( opt.iphys==1 )

  fea = addphys(fea,@navierstokes);     % Add Navier-Stokes equations physics mode.
  fea.phys.ns.eqn.coef{2,end} = { Pr };
  fea.phys.ns.eqn.coef{4,end} = { [num2str(Ra*Pr),'*T'] };
  fea.phys.ns.sfun            = { sf_u sf_u sf_p };   % Set shape functions.

  fea = addphys(fea,@heattransfer);     % Add heat transfer physics mode.
  fea.phys.ht.sfun            = { sf_T };
  fea.phys.ht.eqn.coef{4,end} = { fea.phys.ns.dvar{1} };
  fea.phys.ht.eqn.coef{5,end} = { fea.phys.ns.dvar{2} };
  fea.phys.ht.bdr.sel([ib_l ib_r])  = 1;
  fea.phys.ht.bdr.coef{1,end}{ib_l} = 1;

  fea = parsephys(fea);                 % Check and parse physics modes.

else

  fea.dvar  = { 'u'  'v'  'p'  'T'  };       % Dependent variable name.
  fea.sfun  = { sf_u sf_u sf_p sf_T };       % Shape function.

  % Define equation system.
  cvelx = fea.dvar{1};   % Convection velocities.
  cvely = fea.dvar{2};
  fea.eqn.a.form = { [2 3 2 3;2 3 1 1]      [2;3]                   [1;2] []; ...
                     [3;2]                  [2 3 2 3;2 3 1 1]       [1;3] []; ...
                     [2;1]                  [3;1]                   []    []; ...
                     []                     []                      []    [2 3 2 3;2 3 1 1] };
  fea.eqn.a.coef = { {2*Pr Pr cvelx cvely}   Pr                     -1 []; ...
                     Pr                     {Pr 2*Pr cvelx cvely}   -1 []; ...
                     1                       1                      [] []; ...
                     []                     []                      [] {1 1 cvelx cvely} };
  fea.eqn.f.form = { 1 1 1 1 };
  fea.eqn.f.coef = { 0 [num2str(Ra*Pr),'*',fea.dvar{4}] 0 0 };

  % Define boundary conditions.
  fea.bdr.d = cell(4,n_bdr);
  [fea.bdr.d{1:2,:}] = deal(0);

  fea.bdr.d{4,ib_l}  = 1;
  fea.bdr.d{4,ib_r}  = 0;

  fea.bdr.n = cell(4,n_bdr);
  [fea.bdr.n{4,[ib_t ib_b n_bdr]}] = deal(0);
  [fea.bdr.n{3,setdiff(1:n_bdr,n_bdr)}] = deal(0);
  [fea.bdr.n{1:2,n_bdr}] = deal(0);

end

% Parse and solve problem.
fea       = parseprob(fea);             % Check and parse problem struct.
fea.sol.u = solvestat(fea,'fid',fid);   % Call to stationary solver.

% Postprocessing.
if ( opt.iplot>0 )
  figure
  subplot(1,2,1)
  postplot(fea,'surfexpr','sqrt(u^2+v^2)')
  title('Velocity field')
  subplot(1,2,2)
  postplot(fea,'surfexpr','T')
  title('Temperature')
end

% Error checking.
out.err  = nan;
out.pass = nan;
iref = find( Ra==Ra_ref );
if ( ~isempty(iref) )
  n_evalution_points = 3*l/opt.hmax;
  x_eval = linspace( 0, l, n_evalution_points );
  x_mid  = l/2*ones( 1, n_evalution_points );

  u_eval = evalexpr( 'u', [x_mid; x_eval], fea );
  [u_max,ix] = max( u_eval );
  y_max = x_eval( ix );

  v_eval = evalexpr( 'v', [x_eval; x_mid], fea );
  [v_max,ix] = max( v_eval );
  x_max = x_eval( ix );

  Nu_mean = abs( intbdr( 'Tx', fea, 1, 2 ) + intbdr( 'Tx', fea, 2, 2 ) )/2;

  out.u_max   = u_max;
  out.y_max   = y_max;
  out.v_max   = v_max;
  out.x_max   = x_max;
  out.Nu_mean = Nu_mean;

  out.err = [ abs(u_max_ref(iref)-u_max)/u_max_ref(iref);
              abs(y_max_ref(iref)-y_max)/y_max_ref(iref);
              abs(v_max_ref(iref)-u_max)/v_max_ref(iref);
              abs(x_max_ref(iref)-x_max)/x_max_ref(iref);
              abs(Nu_mean_ref(iref)-Nu_mean)/Nu_mean_ref(iref) ];
  out.pass = all( out.err<opt.tol );
end

if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
% the Variable Span beamformer
% h = sum( gevdVec * gevdVec^H / (mu + lambda)) * PhiX * refVec
% input:    PhiN, (Nch, Nch, Nbin) the noise covariance matrix
%           PhiX, (Nch, Nch, Nbin) the speech covariance matrix
%           mu, the speech distortion/noise reduction trade-off parameter
%           span, 1 <= span <= Nch, dimension of the assumed source space
%           refMic, the reference microphone, default = 1
% output:   h, (Nch, Nbin)  the beamformer coefficients
% % Ziteng Wang @ 201812

function h = VarSpan(PhiX, PhiN, mu, span, refMic)
if nargin < 3
    mu = 1;             %%% typical value {0, 1}
elseif nargin < 4
    span = 1;           %%% default 1 for single-target
elseif nargin < 5
    refMic = 1;
end

[Nch, ~, Nbin] = size(PhiX);
refVec = zeros(Nch, 1);
refVec(refMic) = 1;
h = zeros(Nch, Nbin);

for bin = 1:Nbin
    if rcond(PhiN(:,:,bin)) < eps
        disp(['bin ' num2str(bin) ': Noise covariance ill-conditioned.']);
        PhiN(:,:,bin) = PhiN(:,:,bin) + 1e-10 * eye(Nch);
    end
    [vv, dd] = eig(PhiX(:,:,bin), PhiN(:,:,bin));
    [~, idx] = sort(diag(dd), 'descend');
    tmp = 0;
    for i = 1:span
        tmp = tmp + vv(:,idx(i)) * vv(:,idx(i))' / (mu + dd(idx(i),idx(i)));
    end
    h(:,bin) = tmp * PhiX(:,:,bin) * refVec;
end

author: ZitengWang

repo: MASP

save_path: github-repos/MATLAB/ZitengWang-MASP


text:
function [X,Y,Z] = elec_cm2m(X,Y,Z)

% elec_cm2m - Convert electrode coordinates from centimeters to meters
%
%   [X,Y,Z] = elec_cm2m(X,Y,Z)
%
%   Simply:  X = X / 100;  Y = Y / 100;  Z = Z / 100;
%

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

% Licence:  GNU GPL, no implied or express warranties
% History:  07/2001, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    X = X / 100;  Y = Y / 100;  Z = Z / 100;

author: PatternRecognition

repo: OpenBMI


text:
function matlab_mistake_test01 ( )

%*****************************************************************************80
%
%% MATLAB_MISTAKE_TEST01 puts numbers in an array, except they don't get there!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MATLAB_MISTAKE_TEST01:\n' );
  fprintf ( 1, '  Set X(I) = 1 / ( I - 5 ) for I = 1 to 10\n' );
  fprintf ( 1, '  but set X(5) to zero.\n' );

  x = zeros ( 10, 1 );

  for i = 1 : 10
    if ( x(i) ~= 5 )
      x(i) == 1 / ( i - 5 );
    else
      x(i) = 0;
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now let''s print out X!\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    fprintf ( 1, '  X(%d) = %d\n', i, x(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  That was NOT what we expected!\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


clear; close all; clc;
load '5jia1.mat'
load 'x1'
load 'x2'
load 'x3'
load 'x4'
load 'x5'
load 'uavW.mat'
i = uavW;
clear uavW;

P.gravity = 9.8;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Params for Aersonade UAV
%physical parameters of airframe
P.mass = 13.5;
P.Jx   = 0.8244;
P.Jy   = 1.135;
P.Jz   = 1.759;
P.Jxz  = .1204;
% aerodynamic coefficients
P.S_wing        = 0.55;
P.b             = 2.8956;
P.c             = 0.18994;
P.S_prop        = 0.2027;
P.rho           = 1.2682;
P.k_motor       = 80;
P.k_T_P         = 0;
P.k_Omega       = 0;
P.e             = 0.9;

P.C_L_0         = 0.28;
P.C_L_alpha     = 3.45;
P.C_L_q         = 0.0;
P.C_L_delta_e   = -0.36;
P.C_D_0         = 0.03;
P.C_D_alpha     = 0.30;
P.C_D_p         = 0.0437;
P.C_D_q         = 0.0;
P.C_D_delta_e   = 0.0;
P.C_m_0         = -0.02338;
P.C_m_alpha     = -0.38;
P.C_m_q         = -3.6;
P.C_m_delta_e   = -0.5;
P.C_Y_0         = 0.0;
P.C_Y_beta      = -0.98;
P.C_Y_p         = 0.0;
P.C_Y_r         = 0.0;
P.C_Y_delta_a   = 0.0;
P.C_Y_delta_r   = -0.17;
P.C_ell_0       = 0.0;
P.C_ell_beta    = -0.12;
P.C_ell_p       = -0.26;
P.C_ell_r       = 0.14;
P.C_ell_delta_a = 0.08;
P.C_ell_delta_r = 0.105;
P.C_n_0         = 0.0;
P.C_n_beta      = 0.25;
P.C_n_p         = 0.022;
P.C_n_r         = -0.35;
P.C_n_delta_a   = 0.06;
P.C_n_delta_r   = -0.032;
P.C_prop        = 1.0;
P.M             = 50;
P.epsilon       = 0.1592;
P.alpha0        = 0.4712;

% wind parameters
P.wind_n = 0;
P.wind_e = 0;
P.wind_d = 0;
P.L_u = 200;
P.L_v = 200;
P.L_w = 50;
P.sigma_u = 1.06;
P.sigma_v = 1.06;
P.sigma_w = .7;

% r1 - r8
P.r = P.Jx*P.Jz-P.Jxz^2;
P.r1 = P.Jxz*(P.Jx-P.Jy+P.Jz)/P.r;
P.r2 = P.Jz*(P.Jz-P.Jy)+P.Jxz^2;
P.r3 = P.Jz/P.r;
P.r4 = P.Jxz/P.r;
P.r5 = (P.Jz-P.Jx)/P.Jy;
P.r6 = P.Jxz/P.Jy;
P.r7 = ((P.Jx-P.Jy)*P.Jx+P.Jxz^2)/P.r;
P.r8 = P.Jx/P.r;

% C parameters on p.62
P.C_p_0 = P.r3 * P.C_ell_0 + P.r4 * P.C_n_0;
P.C_p_beta = P.r3 * P.C_ell_beta + P.r4 * P.C_n_beta;
P.C_p_p = P.r3 * P.C_ell_p + P.r4 * P.C_n_p;
P.C_p_r = P.r3 * P.C_ell_r + P.r4 * P.C_n_r;
P.C_p_delta_a = P.r3 * P.C_ell_delta_a + P.r4 * P.C_n_delta_a;
P.C_p_delta_r = P.r3 * P.C_ell_delta_r + P.r4 * P.C_n_delta_r;
P.C_r_0 = P.r4 * P.C_ell_0 + P.r8 * P.C_n_0;
P.C_r_beta = P.r4 * P.C_ell_beta + P.r8 * P.C_n_beta;
P.C_r_p = P.r4 * P.C_ell_p + P.r8 * P.C_n_p;
P.C_r_r = P.r4 * P.C_ell_r + P.r8 * P.C_n_r;
P.C_r_delta_a = P.r4 * P.C_ell_delta_a + P.r8 * P.C_n_delta_a;
P.C_r_delta_r = P.r4 * P.C_ell_delta_r + P.r8 * P.C_n_delta_r;

% compute trim conditions using 'mavsim_chap5_trim.slx'
% initial airspeed
P.Va0 = 35;
gamma = 0; %5*pi/180;  % desired flight path angle (radians)
R     = inf; %150;         % desired radius (m) - use (+) for right handed orbit, can't be 0

% autopilot sample rate
P.Ts = 0.01;
P.tau = 0.5;

% first cut at initial conditions

P.pn0    = 0;  % initial North position
P.pe0    = 0;  % initial East position
P.pd0    = 0;  % initial Down position (negative altitude)

P.u0     = P.Va0; % initial velocity along body x-axis
P.v0     = 0;  % initial velocity along body y-axis
P.w0     = 0;  % initial velocity along body z-axis
P.phi0   = 0;  % initial roll angle
P.theta0 = 0;  % initial pitch angle
P.psi0   = 0;  % initial yaw angle
P.p0     = 0;  % initial body frame roll rate
P.q0     = 0;  % initial body frame pitch rate
P.r0     = 0;  % initial body frame yaw rate

% run trim commands
[x_trim, u_trim]=compute_trim('mavsim_trim',P.Va0,gamma,R);
P.u_trim = u_trim;
P.x_trim = x_trim;

% set initial conditions to trim conditions
% initial conditions
P.pn0    = Xplot2(1,6*i-5);
P.pe0    = Xplot2(1,6*i-4);
P.pd0    = -Xplot2(1,6*i-3);

% P.pn0    = 0;  % initial North position
% P.pe0    = 0;  % initial East position
% P.pd0    = 0;  % initial Down position (negative altitude)

P.u0     = x_trim(4);  % initial velocity along body x-axis
P.v0     = x_trim(5);  % initial velocity along body y-axis
P.w0     = x_trim(6);  % initial velocity along body z-axis
P.phi0   = x_trim(7);  % initial roll angle
P.theta0 = x_trim(8);  % initial pitch angle
P.psi0   = x_trim(9);  % initial yaw angle
P.p0     = x_trim(10);  % initial body frame roll rate
P.q0     = x_trim(11);  % initial body frame pitch rate
P.r0     = x_trim(12);  % initial body frame yaw rate

% compute different transfer functions
[T_phi_delta_a,T_chi_phi,T_theta_delta_e,T_h_theta,T_h_Va,T_Va_delta_t,T_Va_theta,T_v_delta_r]...
    = compute_tf_model(x_trim,u_trim,P);

% linearize the equations of motion around trim conditions
[A_lon, B_lon, A_lat, B_lat] = compute_ss_model('mavsim_trim',x_trim,u_trim);

eig_lon = eig(A_lon);
eig_lat = eig(A_lat);

compute_gains;

% longitudinal state-machine parameters
P.altitude_take_off_zone = 50;
P.altitude_hold_zone = 10;

% sensor parameters(chapter 7)
P.sigma_gyro = 0.13;
P.sigma_accel = 0.0025;
P.beta_abs_pres = 0.125;
P.sigma_abs_pres = 0.01;
P.beta_diff_pres = 0.02;
P.sigma_diff_pres = 0.002;

P.Ts_gps = 1;
P.k_gps = 1/1000;
P.sigma_gps_n = 0.21;
P.sigma_gps_e = 0.21;
P.sigma_gps_altitude = 0.4;
P.sigma_gps_V_g = 0.05;

% filter parameters(chapter 8)
P.alpha_lpf_gyro = 0.1;
P.alpha_lpf_static_pres = 0.1;
P.alpha_lpf_diff_pres = 0.1;

P.bias_gyro_x = 0;             % x-gyro bias
P.bias_gyro_y = 0;             % y-gyro bias
P.bias_gyro_z = 0;             % z-gyro bias

% autopilot guidance model coefficients(chapter 9)
wn_h = 1.5;
zeta_h = 0.707;
P.b_hdot = 2 * zeta_h * wn_h;
P.b_h = wn_h^2;

wn_chi = 0.7;
zeta_chi = 0.707;
P.b_chidot = 2 * zeta_chi * wn_chi;
P.b_chi = wn_chi^2;

P.b_Va = 10;
P.gamma_max = 45*pi/180;

% chapter 10 parameters
P.k_path = 0.05;
P.chi_inf = 40*pi/180;
P.k_orbit = 0.01;

% chapter 11 - path manager
% number of waypoints in data structure
P.size_waypoint_array = 100;
P.R_min = P.Va0^2/P.gravity/tan(P.roll_max);

% create random city map
city_width      = 2000;  % the city is of size (width)x(width)
building_height = 300;   % maximum height of buildings
%building_height = 1;   % maximum height of buildings (for camera)
num_blocks      = 5;    % number of blocks in city
street_width    = .8;   % percent of block that is street.
P.h0 = 100;
P.pd0           = -P.h0;  % initial height of MAV
P.map = createWorld(city_width, building_height, num_blocks, street_width);


text:
function [a,r,f] = femq1_diff(xy,ev)
%femq1_diff  vectorized bilinear coefficient matrix generator
%   [A,Q,f] = femq1_diff(xy,ev);
%   input
%          xy         vertex coordinate vector
%          ev         element mapping matrix
%   output
%          A          stiffness matrix
%          Q          mass matrix
%          f          rhs vector
%
%   Natural boundary conditions apply. Dirichlet conditions
%   must be explicitly enforced by calling function nonzerobc.
%   IFISS function: DJS; 4 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
x=xy(:,1); y=xy(:,2);
nvtx=length(x);
nel=length(ev(:,1));
lx=max(x)-min(x); ly=max(y)-min(y);
hx=max(diff(x)); hy=max(diff(y));
fprintf('setting up Q1 diffusion matrices...  ')
%
% initialise global matrices
      a = sparse(nvtx,nvtx);
      r = sparse(nvtx,nvtx);
      f = zeros(nvtx,1);
%
% set up 2x2 Gauss points
      gpt=1.0e0/sqrt(3.0e0);
      s(1) = -gpt;  t(1) = -gpt;
      s(2) =  gpt;  t(2) = -gpt;
      s(3) =  gpt;  t(3) =  gpt;
      s(4) = -gpt;  t(4) =  gpt;
%
% inner loop over elements
        for ivtx = 1:4
        xl_v(:,ivtx) = x(ev(:,ivtx));
        yl_v(:,ivtx) = y(ev(:,ivtx));
	    end
        ae = zeros(nel,4,4);
        re = zeros(nel,4,4);
        fe = zeros(nel,4);
%  loop over 2x2 Gauss points
         for igpt = 1:4
         sigpt=s(igpt);
         tigpt=t(igpt);
%  evaluate derivatives etc
         [jac,invjac,phi,dphidx,dphidy] = deriv(sigpt,tigpt,xl_v,yl_v);
         rhs = gauss_source(sigpt,tigpt,xl_v,yl_v);
            for j = 1:4
               for i = 1:4
               ae(:,i,j) = ae(:,i,j)  + dphidx(:,i).*dphidx(:,j) .* invjac(:);
               ae(:,i,j) = ae(:,i,j)  + dphidy(:,i).*dphidy(:,j) .* invjac(:);
               re(:,i,j) = re(:,i,j)  + phi(:,i).*phi(:,j) .* jac(:);
               end
            fe(:,j) = fe(:,j)  + rhs(:) .* phi(:,j) .* jac(:);
	        end
% end of Gauss point loop
         end
% perform assembly of global matrix  and source vector
      for krow=1:4
	  nrow=ev(:,krow);
          for kcol=1:4
		  ncol=ev(:,kcol);
          a = a + sparse(nrow,ncol,ae(:,krow,kcol),nvtx,nvtx);
          r = r + sparse(nrow,ncol,re(:,krow,kcol),nvtx,nvtx);
          end
      f(nrow,1) = f(nrow,1) + fe(:,krow);
      end
%
fprintf('done\n')
return


text:
%%%%%%%VERSION 3
%%ANOTHER DESCRIBTION OF GABOR FILTER

%The Gabor filter is basically a Gaussian (with variances sx and sy along x and y-axes respectively)
%modulated by a complex sinusoid (with centre frequencies U and V along x and y-axes respectively)
%described by the following equation
%%
%               1                -1     x  ^    y  ^
%%% Gi(x,y) = ---------- * exp ([----{(----) 2+(----) 2}])*Mi(x,y,f);
%            2*pi*sx*sy           2    sx       sy
%%% i =1,2
%%% M1(x,y,f) = cos[2*pi*f*sqrt(x^2+y^2)];
%%% M2(x,y,f) = cos[2*pi*f*(x*cos(theta) + y*sin(theta)];

%% Describtion :

%% I : Input image
%% Sx & Sy : Variances along x and y-axes respectively
%% f : The frequency of the sinusoidal function
%% theta : The orientation of Gabor filter

%% G1 & G2 : The output filters as described above
%% gabout1 & gabout2 : The output filtered images

%%  Author : Ahmad poursaberi  e-mail : a.poursaberi@ece.ut.ac.ir
%%          Faulty of Engineering, Electrical&Computer Department,Tehran
%%          University,Iran,June 2004

function [G1,G2,gabout1,gabout2] = gaborfilter(I,Sx,Sy,f,theta);

if isa(I,'double')~=1
    I = double(I);
end

for x = -fix(Sx):fix(Sx)
    for y = -fix(Sy):fix(Sy)
        M1 = cos(2*pi*f*sqrt(x^2+y^2));
        M2 = cos(2*pi*f*(x*cos(theta)+y*sin(theta)));
        G1(fix(Sx)+x+1,fix(Sy)+y+1) = (1/(2*pi*Sx*Sy)) * exp(-.5*((x/Sx)^2+(y/Sy)^2))*M1;
        G2(fix(Sx)+x+1,fix(Sy)+y+1) = (1/(2*pi*Sx*Sy)) * exp(-.5*((x/Sx)^2+(y/Sy)^2))*M2;
    end
end

Imgabout1 = conv2(I,double(imag(G1)),'same');
Regabout1 = conv2(I,double(real(G1)),'same');

Imgabout2 = conv2(I,double(imag(G2)),'same');
Regabout2 = conv2(I,double(real(G2)),'same');

gabout1 = sqrt(Imgabout1.*Imgabout1 + Regabout1.*Regabout1);
gabout2 = sqrt(Imgabout2.*Imgabout2 + Regabout2.*Regabout2);

author: Sable

repo: mcbench-benchmarks


text:
function output = calc_traversal_dist(ai)

% This function will generate position coordinates of chain code (ai). Number of
% harmonic elements (n), and number of points for reconstruction (m) must be
% specified.

    x_ = 0;
    y_ = 0;

    for i = 1 : size(ai, 2)
        x_ = x_ + sign(6 - ai(i)) * sign(2 - ai(i));
        y_ = y_ + sign(4 - ai(i)) * sign(ai(i));
        p(i, 1) = x_;
        p(i, 2) = y_;
    end

    output = p;

end


text:
function [yd3] = oz2yd3(oz)
% Convert volume from US liquid ounces to cubic yards.
% Chad Greene 2012
yd3 = oz*0.000038680716307;


text:
function [check,maxerr,yfit] = cvxfitfeas(U,y,err)

% Checks if a data set can be approximated with a convex function subject
% to a maximum residual constraint.
%
% Notation:
%
% U - matrix of input data points, with rows denoting the different points
% y - outputs corresponding to the inputs (row vector)
% err - maximum allowable residual
% check - returns 1 if it is possible to fit function, 0 otherwise
% maxerr - maximum residual of convex fit
% yfit - output values defining the fit

m = length(y);
n = length(U(1,:));

A1 = [-eye(m) zeros(m,m*n) -ones(m,1);
    eye(m) zeros(m,m*n) -ones(m,1)];
b1 = [-y'; y'];

A2 = zeros(m*(m-1),m+m*n+1);
lc = 1;
for i = 1:m
    for j = 1:m
        if i ~= j
            A2(lc,i) = 1;
            A2(lc,j) = -1;
            A2(lc,[1+m+(i-1)*n:m+i*n]) = U(j,:)-U(i,:);
            lc = lc + 1;
        end
    end
end

b2 = zeros(m*(m-1),1);

A = [A1;A2];
b = [b1;b2];

[xopt,maxerr,exit] = linprog([zeros(1,m+m*n) 1],A,b,[],[],[],[],[],optimset('disp','iter'));
if exit <= 0
    disp('Warning - Optimization did not terminate property.');
end
yfit = xopt(1:m);

if maxerr <= err
    check = 1;
else
    check = 0;
end

end


text:
% op_rmNworstaverages.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [out,metric,badAverages]=op_rmNworstaverages(in,n);
%
% DESCRIPTION:
% Removes motion corrupted averages from a dataset containing multiple
% averages.  The N most badly motion corrupted averages are discarded.
%
% INPUTS:
% in         = input data in matlab structure format
% n          = number of bad averages to remove
%
% OUTPUTS:
% out         = Output dataset following removal of motion corrupted averages.
% metric      = Vector of unlikeness metrics corresponding to all input
%               averages.
% badAverages = Indices of the averages that were removed.

function [out,metric,badAverages]=op_rmNworstaverages(in,n);

if in.flags.averaged
    error('ERROR:  Averaging has already been performed!  Aborting!');
end

if ~in.flags.addedrcvrs
    error('ERROR:  Receivers should be combined first!  Aborting!');
end

%first, make a metric by subtracting all averages from the first average,
%and then taking the sum of all all the spectral points.
if in.dims.subSpecs>0
    SS=in.sz(in.dims.subSpecs);
else
    SS=1;
end
infilt=op_filter(in,10);
%inavg=op_averaging(infilt);
inavg=op_median(infilt);
for k=1:in.sz(in.dims.averages)
    for m=1:SS
            metric(k,m)=sum((real(infilt.specs(:,k,m))-(real(inavg.specs(:,m)))).^2);
    end
end

%find the average and standard deviation of the metric
avg=mean(metric);
stdev=std(metric);

%Now z-transform the metric so that it is centered about zero, and they
%have a standard deviation of 1.0.
zmetric=(metric-avg)/stdev;

for m=1:SS
    P(m,:)=polyfit([1:in.sz(in.dims.averages)]',zmetric(:,m),2);
    figure
    plot([1:in.sz(in.dims.averages)],zmetric(:,m),'.',[1:in.sz(in.dims.averages)],polyval(P(m,:),[1:in.sz(in.dims.averages)]));
end

%Now make a mask that represents the locations of the averages
%whose metric values are more than nsd standard deviations away from the
%mean metric value.

%first sort the zmetric array to find the n highest values:

[zmetric_sorted,inds]=sort(zmetric-polyval(P,[1:in.sz(in.dims.averages)])',1,'descend');

mask=zeros(size(zmetric));

for l=1:SS
    %mask(:,l)=metric(:,l)>(avg(l)+(nsd*stdev(l))) | metric(:,l)<(avg(l)-(nsd*stdev(l)));
    %mask(:,l)=metric(:,l)>(polyval(P(l,:),[1:in.sz(in.dims.averages)])'+(nsd*stdev(l))) | metric(:,l)<(polyval(P(l,:),[1:in.sz(in.dims.averages)])'-(nsd*stdev(l)));
    %mask(:,l)=metric(:,l)>(polyval(P(l,:),[1:in.sz(in.dims.averages)])'+(nsd*stdev(l)));
    %mask(:,l)=(metric(:,l)-polyval(P(l,:),[1:in.sz(in.dims.averages)])')==max((metric(:,l)-polyval(P(l,:),[1:in.sz(in.dims.averages)])'));
    for b=1:n
        mask(inds(b),l)=1;
    end
end

%Unfortunately, if one average is corrupted, then all of the subspecs
%corresponding to that average have to be thrown away.  Therefore, take the
%minimum intensity projection along the subspecs dimension to find out
%which averages contain at least one corrupted subspec:
if size(mask,2)>1
    mask=sum(mask')'>0;
end

%now the corrupted and uncorrupted average numbers are given by:
badAverages=find(mask);
goodAverages=find(~mask);

%make a new fids array containing only good averages
fids=in.fids(:,goodAverages,:,:);

%%re-calculate Specs using fft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%re-calculate the sz variable
sz=size(fids);

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.averages=length(goodAverages) * in.rawSubspecs;

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
%%*******************************************************************
%% Prod2: compute the block diagonal matrix A*B
%%
%%    C = Prod2(blk,A,B,options);
%%
%% INPUT:  blk = a cell array describing the block structure of A and B
%%         A,B = square matrices or column vectors.
%%
%%   options = 0  if no special structure
%%             1  if C is symmetric
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function  C = Prod2(blk,A,B,options);

   global spdensity

   if (nargin == 3); options = 0; end;
   iscellA = iscell(A); iscellB = iscell(B);
%%
   if (~iscellA & ~iscellB)
      if (size(blk,1) > 1);
         error('Prod2: blk and A,B are not compatible');
      end;
      if strcmp(blk{1},'s')
         numblk = length(blk{2});
         isspA = issparse(A); isspB = issparse(B);
         if (numblk > 1)
            if ~isspA; A=sparse(A); isspA=1; end
            if ~isspB; B=sparse(B); isspB=1; end
         end
         %%use_matlab = (options==0 & ~isspA & ~isspB) | (isspA & isspB);
         use_matlab = (~isspA & ~isspB) | (isspA & isspB);
         if (use_matlab)
            C = A*B;
            if (options==1); C = 0.5*(C+C'); end;
         else
            C = mexProd2(blk,A,B,options);
         end
         checksparse = (numblk==1) & (isspA | isspB);
         if (checksparse)
            n2 = sum(blk{2}.*blk{2});
            if (mexnnz(C) <= spdensity*n2);
               if ~issparse(C); C = sparse(C); end;
            else
               if issparse(C); C = full(C); end;
            end
         end
      elseif (strcmp(blk{1},'q') | strcmp(blk{1},'l') | strcmp(blk{1},'u'))
         C = A.*B;
      end
   else
      error(['Prod2: A,B must be matrices']);
   end
%%*******************************************************************

author: yu-jiang

repo: radpbook

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


text:
function yth = SB_CoarseGrain(y,howtocg,numGroups)
% SB_CoarseGrain   Coarse-grains a continuous time series to a discrete alphabet.
%
%---INPUTS:
% howtocg, the method of coarse-graining
%
% numGroups, either specifies the size of the alphabet for 'quantile' and 'updown'
%       or sets the time delay for the embedding subroutines

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
% Check inputs, preliminaries:
% ------------------------------------------------------------------------------

% Quantile puts an equal number into each bin
if nargin < 3
    howtocg = 'quantile';
end

N = length(y); % length of the input sequence

if ~ismember(howtocg,{'updown','quantile','embed2quadrants','embed2octants'})
    error('Unknown coarse-graining method ''%s''',howtocg);
end

% ------------------------------------------------------------------------------
% Some coarse-graining/symbolization methods require initial processing:
% ------------------------------------------------------------------------------
switch howtocg
case 'updown'
   y = diff(y);
   N = N - 1; % the time series is one value shorter than the input because of differencing
   howtocg = 'quantile'; % successive differences and then quantiles

case {'embed2quadrants','embed2octants'}
	% Construct the embedding

	if strcmp(numGroups,'tau')
        % First zero-crossing of the autocorrelation function
        tau = CO_FirstCrossing(y,'ac',0,'discrete');
    else
        tau = numGroups;
	end
	if tau > N/25; tau = floor(N/25); end
	m1 = y(1:end-tau);
	m2 = y(1+tau:end);

	% Look at which points are in which angular 'quadrant'
	upr = find(m2 >= 0); % points above the axis
	downr = find(m2 < 0); % points below the axis

	q1r = upr(m1(upr) >= 0); % points in quadrant 1
	q2r = upr(m1(upr) < 0); % points in quadrant 2
	q3r = downr(m1(downr) < 0); % points in quadrant 3
	q4r = downr(m1(downr) >= 0); % points in quadrant 4
end

% ------------------------------------------------------------------------------
% Do the coarse graining
% ------------------------------------------------------------------------------
switch howtocg
    case 'quantile'
        th = quantile(y,linspace(0,1,numGroups+1)); % thresholds for dividing the time-series values
        th(1) = th(1)-1; % this ensures the first point is included
        % turn the time series into a set of numbers from 1:numGroups
        yth = zeros(N,1);
        for i = 1:numGroups
            yth(y > th(i) & y <= th(i+1)) = i;
        end

    case 'embed2quadrants' % divides based on quadrants in a 2-D embedding space
		% create alphabet in quadrants -- {1,2,3,4}
		yth = zeros(length(m1),1);
		yth(q1r) = 1; yth(q2r) = 2; yth(q3r) = 3; yth(q4r) = 4;

	case 'embed2octants' % divide based on octants in 2-D embedding space
		o1r = q1r(m2(q1r)<m1(q1r)); % points in octant 1
		o2r = q1r(m2(q1r)>=m1(q1r)); % points in octant 2
		o3r = q2r(m2(q2r)>=-m1(q2r)); % points in octant 3
		o4r = q2r(m2(q2r)<-m1(q2r)); % points in octant 4
		o5r = q3r(m2(q3r)>=m1(q3r)); % points in octant 5
		o6r = q3r(m2(q3r)<m1(q3r)); % points in octant 6
		o7r = q4r(m2(q4r)<-m1(q4r)); % points in octant 7
		o8r = q4r(m2(q4r)>=-m1(q4r)); % points in octant 8

		% create alphabet in octants -- {1,2,3,4,5,6,7,8}
		yth = zeros(length(m1),1);
		yth(o1r) = 1; yth(o2r) = 2; yth(o3r) = 3; yth(o4r) = 4;
		yth(o5r) = 5; yth(o6r) = 6; yth(o7r) = 7; yth(o8r) = 8;
end

if any(yth == 0)
    error('All values in the sequence were not assigned to a group')
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%  Figure 10.18      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  fig10_18.m is a script to generate Fig 10.18,
%  the symmetric  rootlocus for the LQR symmetric rootlocus
%  compensator of the satellite position control, non-collocated case

% parameter values
m=[1, 0.1]; k=[0, 0.091] ; d=[0, 0.0036]; k1=[0, 0.4];
% call model
[f,g,h,j] = twomass(m,k,d);
s=[f, g;h, 0];r=[0*g;1]; n=s\r;nx=n(1:4);nu=n(5);
% call model
[f1,g,h,j] = twomass(m,k1,d);
% Form G(s)G(-s) in state-space
a=[f, 0*f;
-h'*h, -f'];
b=[g;0*g];
c=[0*h, g'];
d=[0];

hold off; clf;
sysGG=ss(a,b,c,d);
rlocus(sysGG);
grid;
v=[-2 2 -1.5 1.5];axis(v);
title('Fig. 10.18 Symmetric rootlocus for the satellite')


text:
function [AlgOption,TSP] = InterfaceMMAS(TSPfile,AntNum,alpha,beta,rho,MaxITime)
%:Get the input parameters for MMAS
%Read in tsp data
[Dimension,Nodes,Weights,Name] = GetTSPData(TSPfile);
fprintf('Data of problem-%s have been read in!\n',Name);
%Set TSP datas
TSP = InitProblem(Dimension,Nodes,Weights,Name);
%Set parameters for MMAS algorithm
AlgOption = InitParameter(Dimension,AntNum,alpha,beta,rho,MaxITime);
% AlgOption = InitParameter(Dimension,Dimension,alpha,beta,rho,MaxITime);

%% --------------------------------------------------------------
function AlgorithmParas = InitParameter(Dimension,AntNum,alpha,beta,rho,MaxITime)
AlgorithmParas.n = Dimension; % nodes number in TSP
AlgorithmParas.m = AntNum; % ants number
AlgorithmParas.alpha = alpha; % pheromeno exponential
AlgorithmParas.beta = beta; % heuristic exponential
AlgorithmParas.rho = rho; % vapor parameter
AlgorithmParas.MaxITime = MaxITime; % Maximum Iterative Time
AlgorithmParas.delta = 0.05; % Parameter of Pheromone Trail Smoothing(PTS): (0,1) Turn Off By Setting 0
AlgorithmParas.lambda = 0; % coefficient of Average Node Branching: (0,1) Turn Off By Setting 0
AlgorithmParas.ANBmin = 2; % coefficient of minimum ANB
AlgorithmParas.DispInterval = 5; % Display Interval: Turn Off Display By Setting 0
AlgorithmParas.ReStartCount = 50; % Restart while no improvment for 50 iterations
%% --------------------------------------------------------------
function Problem = InitProblem(Dimension,Nodes,NodeWeight,Name)
n = Dimension;
MatrixTau = ones(n,n)-eye(n,n);
Problem.nodes = Nodes;
Problem.weights = NodeWeight;
Problem.tau = MatrixTau;
Problem.name = Name;


text:
% eml_em_test.m
% compare aspire and matlab E-ML-EM
% Copyright Jan 1998	Jeff Fessler, The University of Michigan

%
% generate data
%
if ~isvar('yi'), printm 'setup for eml_em_test'
	if has_aspire
		f.dir	= test_dir;
		f.wtf	= [f.dir 't,g.wtf'];
		f.yi	= [f.dir 'yi.fld'];
		f.ci	= [f.dir 'ci.fld'];
		f.ri	= [f.dir 'ri.fld'];
		f.mask	= [f.dir 'mask.fld'];
	end
	em_test_setup
prompt
end

%
% matlab iterations
%
if ~isvar('xmat'), printm 'matlab E-ML-EM'
	f.niter = 9;
	xinit = ig.ones; % uniform
	xmat = eml_em(xinit(ig.mask), G, yi(:), ci(:), ri(:), [], f.niter);
	xmat = ig.embed(xmat);

	im clf, im(xmat, 'matlab E-ML-EM iterates')
prompt
end

if ~has_aspire, return, end

%
% aspire iterations - this is for Fessler's testing only!
%
if ~isvar('xasp'), printm 'aspire E-ML-EM'

	f.init	= [f.dir 't,init.fld'];
	f.out = [f.dir 't,out.fld'];
	fld_write(f.init, xinit, 'check', 0)
	if (exist(f.out) == 2), delete(f.out), end

	f.alg = 'em,1';
	f.saver = '-';
	f.saver = 'stack,1';
	f.penal = '-';
	f.method = sprintf('@%d@%s@%s', f.niter-1, f.alg, f.penal);
	f.scaleinit = 0;

	if 1
		f.com = sprintf('i -chat 0 empl2 %s %s  %s %s %s 1 %s %s  %s %s 0 1e9 %d -', ...
			f.out, f.init, f.yi, f.ci, f.ri, f.wtf, ...
			f.mask, f.method, f.saver, f.scaleinit);
		os_run(f.com)

		xasp = double(fld_read(f.out));
	end
end

if 1
	im(221, xmat, 'xhat matlab')
	im(222, xasp, 'xhat aspire')
	im(223, (xasp-xmat)/max(xmat(:)), 'aspire-matlab')

	t = vcorrcoef(xasp, xmat);
	printf('corr. %g,%g', t, t-1)

	t1 = eql_obj(xmat, G, yi(:), ci(:), ri(:), [], ig.mask);
	t2 = eql_obj(xasp, G, yi(:), ci(:), ri(:), [], ig.mask);

	if im
		subplot(224)
		plot(0:f.niter-1, t1-t1(1), '-o', 0:f.niter-1, t2-t1(1), '-x')
		legend('mat', 'asp', 4), xlabel iteration, ylabel objective
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ n_data, x, fx ] = tran06_values ( n_data )

%*****************************************************************************80
%
%% TRAN06_VALUES returns some values of the order 6 transportation function.
%
%  Discussion:
%
%    The function is defined by:
%
%      TRAN06(x) = Integral ( 0 <= t <= x ) t^6 * exp(t) / ( exp(t) - 1 )^2 dt
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%    special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
     0.56843405953641209574E-14, ...
     0.59601180165247401484E-08, ...
     0.60978424397580572815E-05, ...
     0.61578909866319494394E-02, ...
     0.18854360275680840514E+00, ...
     0.13319251347921659134E+01, ...
     0.50857202271697616755E+01, ...
     0.13729222365466557122E+02, ...
     0.29579592481641441292E+02, ...
     0.88600835706899853768E+02, ...
     0.10916037113373004909E+03, ...
     0.18224323749575359518E+03, ...
     0.23765383125586756031E+03, ...
     0.29543246745959381136E+03, ...
     0.50681244381280455592E+03, ...
     0.63878231134946125623E+03, ...
     0.72699203556994876111E+03, ...
     0.73230331643146851717E+03, ...
     0.73248692015882096369E+03, ...
     0.73248700462879996604E+03 ];

  x_vec = [ ...
       0.0019531250E+00, ...
       0.0312500000E+00, ...
       0.1250000000E+00, ...
       0.5000000000E+00, ...
       1.0000000000E+00, ...
       1.5000000000E+00, ...
       2.0000000000E+00, ...
       2.5000000000E+00, ...
       3.0000000000E+00, ...
       4.0000000000E+00, ...
       4.2500000000E+00, ...
       5.0000000000E+00, ...
       5.5000000000E+00, ...
       6.0000000000E+00, ...
       8.0000000000E+00, ...
      10.0000000000E+00, ...
      15.0000000000E+00, ...
      20.0000000000E+00, ...
      30.0000000000E+00, ...
      50.0000000000E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function Y = svmlfwd(net, X, Y)
% SVMLFWD - Wrapper for SVMlight: Prediction
%
%   YPRED = SVMLFWD(NET, X)
%   Predict labels YPRED for points X using the SVMlight stored in NET.
%   Each row of X is one data point, X is of size [N D] for N points of
%   dimensionality D. The result Y is a column vector of size [N 1]. NET
%   is a structure holding information on the trained SVMlight model, as
%   it is output by SVMLTRAIN. In particular, the SVMlight model stored
%   in file NET.fname is used to classify the points X.
%   YPRED = SVMLFWD(NET, X, Y) also provides the actual target for the
%   test data. When using this syntax, SVMlight will print out the
%   accuracy on the test set, as well as precision and recall.
%   SVMLFWD(NET, FNAME), where FNAME is the name of an SVMlight data
%   file, computes the predictions for test data that are stored in file
%   FNAME.
%
%   See also SVML, SVMLTRAIN, SVM_CLASSIFY
%

%
% Copyright (c) by Anton Schwaighofer (2002)
% $Revision: 1.2 $ $Date: 2002/02/19 12:23:16 $
% mailto:anton.schwaighofer@gmx.net
%
% This program is released unter the GNU General Public License.
%

error(nargchk(2, 3, nargin));
error(consist(net, 'svml'));

if nargin<3,
  Y = [];
end
fname = net.fname;
if ischar(X),
  testdata = X;
  deleteData = 0;
else
  testdata = [fname '.testdata'];
  svmlwrite(testdata, X, Y);
  deleteData = 1;
end
preddata = [fname '.preddata'];
status = svm_classify(net.options, testdata, net.fname, preddata);
Y = svmlread(preddata);
delete(preddata);
if deleteData,
  delete(testdata);
end
if status~=0,
  error(sprintf('Error when calling SVMlight. Status = %i', status));
end

author: viggin

repo: domain-adaptation-toolbox


text:
function [ln,lf]=dfa(varargin)
%
% [ln,lf]=dfa(x,p,integrateFlag,minBoxSize,maxBoxSize,slideWindowFlag)
%
%
% Wrapper to the DFA Algorithm in:
%    http://www.physionet.org/physiotools/dfa/
%
% References:
% Peng C-K, Buldyrev SV, Havlin S, Simons M, Stanley HE, Goldberger AL. Mosaic organization of DNA nucleotides. Phys Rev E 1994;49:1685-1689.
% Peng C-K, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos 1995;5:82-87.
%
% Please cite at least one of the above publications when referencing this
% material.
%
% Required Input Options are:
%
% x
%       A Nx1 vector of doubles.
%
% Optional Input Options are:
%
% p
%       Detrend using a polynomial of degree p (default: p=1, linear
%       detrending).
%
% integrateFlag
%       Input series is already integrated ( default= false ).
%
% minBoxSize
%       Smallest box width (default: 2p+2)
%
% maxBoxSize
%       Largest box width (default: N/4)
%
% slideWindowFlag
%       Sliding window DFA (default =false);
%
%
% The Output variables are:
%
% ln
%       A (MaxBoxSize -MinBoxSize) x 1 vector of log(boxsize)
%
% lf
%       A (MaxBoxSize -MinBoxSize) x 1 vector of the log of the root
%       mean square fluctuation for the given boxsize.
%
%
% Written by Ikaro Silva, 2014
% Last Modified: November 21, 2014
% Version 1.0
%
% Since 0.9.8
%
% %Example:
%
%  gqrs('mitdb/117');
%  [rr]=ann2rr('mitdb/117','qrs');
%  [ln,lf]=dfa(rr);
%  plot(ln,lf)

%
%
% See also MSENTROPY, SURROGATE

%endOfHelp

persistent javaWfdbExec config
if(isempty(javaWfdbExec))
    [javaWfdbExec,config]=getWfdbClass('dfa');
end

%Set default pararamter values
inputs={'x','p','integrateFlag','minBoxSize','maxBoxSize','slideWindowFlag'};
p=[];
integrateFlag=false;
minBoxSize=[];
maxBoxSize=[];
slideWindowFlag=false;
wfdb_argument={};
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};'])
    end
end
if(~isempty(p))
    wfdb_argument{end+1}='-d';
    wfdb_argument{end+1}=[num2str(p)];
end
if(integrateFlag)
    wfdb_argument{end+1}='-i';
end
if(~isempty(minBoxSize))
    wfdb_argument{end+1}='-l';
    wfdb_argument{end+1}=[num2str(minBoxSize)];
end
if(~isempty(maxBoxSize))
    wfdb_argument{end+1}='-u';
    wfdb_argument{end+1}=[num2str(maxBoxSize)];
end
if(slideWindowFlag)
    wfdb_argument{end+1}='-s';
end

javaWfdbExec.setArguments(wfdb_argument);

if(config.inOctave)
    x=cellstr(num2str(x));
    x=java2mat(javaWfdbExec.execWithStandardInput(x));
    Nx=x.size;
    out=cell(Nx,1);
    for n=1:Nx
        out{n}=x.get(n-1);
    end
else
    out=cell(javaWfdbExec.execWithStandardInput(x).toArray);
end

M=length(out);
ln=zeros(M,1)+NaN;
lf=zeros(M,1)+NaN;
for m=1:M
    str=out{m};
    sep=regexp(str,'\s');
    ln(m)=str2num(str(1:sep));
    lf(m)=str2num(str(sep(1):end));
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function tests = test_spm_dcm_bmr_all
% Unit Tests for test_spm_dcm_bmr_all
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_dcm_bmr_all.m 6915 2016-11-01 17:38:08Z peter $

tests = functiontests(localfunctions);

% -------------------------------------------------------------------------
function test_bmr_all(testCase)

data_path = get_data_path();

% Load PEB of full DCM
PEB = load(fullfile(data_path,'PEB_test.mat'));
PEB = PEB.PEB;
PEB = PEB(1);

% Load DCMs
GCM = load(fullfile(data_path,'models','GCM_simulated.mat'));
GCM = GCM.GCM;

% Prune connections on the B-matrix only
[rPEB,BMR,BMA] = spm_dcm_bmr_all(PEB,{'B'});

% Check the only surviving B-matrix (group difference) parameter is R1->R2
b_reduced = BMA.Cp(:,2);
expected  = logical([0 0 0 0 1 0]);
testCase.assertTrue(all(b_reduced(~expected) < 1e-4));
testCase.assertTrue(all(b_reduced(expected) > 1e-4));

% -------------------------------------------------------------------------
function data_path = get_data_path()

data_path = fullfile( spm('Dir'), 'tests', ...
    'data', 'fMRI', 'simulated_2region');

author: spm


text:
% Make the following HHMM
%
%     LH                  RH
%    /                      \
%   /                        \
%  LR -> UD -> RL -> DU       RL -> UD -> LR -> DU
%   \
%    \
%     Q1 -> Q2
%
% where level 1 is fully interconnected (not shown)
% level 2 is left-right
% and each model at level 3 is a 2 state LR shared HMM

Qsizes = [2 4 2];
D = 3;

% LEVEL 1

startprob1 = 'ergodic';
transprob1 = 'ergodic';

% LEVEL 2

startprob = zeros(2, 4);
%        Q1  Q2
startprob(1, 1) = 1;
startprob(2, 3) = 1;

transprob = zeros(2, 4, 4);
transprob(1,:,:) = [0 1 0 0
		    0 0 1 0
		    0 0 0 1
		    0 0 0 1];
transprob(2,:,:) = [0 0 0 1
		    1 0 0 0
		    0 1 0 0
		    0 0 0 1];

Q2args = {'startprob', startprob, 'transprob', transprob};

% always terminate in state 4 (default)
% F2args

% LEVEL 3

% Defaults are fine: always start in state 1, left-right model, finish in state 2

% OBS LEVEl

chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
Osize = length(chars);

obsprob = zeros([4 2 Osize]);
%       Q2 Q3 O
obsprob(1, 1, find(chars == 'L')) =  1.0;
obsprob(1, 2, find(chars == 'l')) =  1.0;

obsprob(2, 1, find(chars == 'U')) =  1.0;
obsprob(2, 2, find(chars == 'u')) =  1.0;

obsprob(3, 1, find(chars == 'R')) =  1.0;
obsprob(3, 2, find(chars == 'r')) =  1.0;

obsprob(4, 1, find(chars == 'D')) =  1.0;
obsprob(4, 2, find(chars == 'd')) =  1.0;

Oargs = {'CPT', obsprob};

bnet = mk_hhmm3('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', 1, 'Oargs', Oargs, 'Q1args', Q1args, 'Q2args', Q2args);

T = 20;
usecell = 0;
evidence = sample_dbn(bnet, T, usecell);
%chars(evidence(end,:))

Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6;
Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];

pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, obs, chars);

eclass = bnet.equiv_class;
S=struct(bnet.CPD{eclass(Q2,2)})

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [Nbegcol,Nlencol,Nrowndx] = lindosparse(cJacobian);

Nbegcol = [];
Nrowndx = [];
Nlencol = [];
top = 0;
for i = 1:size(cJacobian,2)
    [ii,jj,kk] = find(cJacobian(:,i));
    if isempty(ii)
        Nbegcol = [Nbegcol top];
        Nlencol = [Nlencol 0];
    else
        Nbegcol = [Nbegcol top];
        Nrowndx = [Nrowndx ii(:)'-1];
        Nlencol = [Nlencol length(ii)];
        top = top + length(ii);
    end
end
if  isempty(Nrowndx)
    Nrowndx = [];
end
Nbegcol = [Nbegcol sum(Nlencol)];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = iff(varargin)
%IFF Logical equivalence
%
% IFF(X,Y) creates a mixed integer representation of
% the constraint X <--> Y, i.e. Y is true iff X is true.
%
% Syntax
%   F = iff(X,Y)
%
% Input
%   X : binary SDPVAR variable or a set of linear (in)equalities
%   Y : binary SDPVAR variable or a set of linear (in)equalities
%
% Output
%   F : Constraint object
%
% Examples
%
%  binvar X,Y; F = iff(X,Y);
%  sdpvar X;binvar Y; F = iff(X>=5,Y);
%  sdpvar X;binvar Y; F = iff(Y,X==5);
%
% Overloading
%
% The iff overloads == for logic constraints.
%
%  sdpvar X;binvar Y; F = ((X>=5) == Y);
%  sdpvar X;binvar Y; F = (Y == (X==5));
%
%
% Note
%  The function IFF is not complete, but will be
%  improved upon in future releases.
%
%   See also @SDPVAR/AND, @SDPVAR/OR, IMPLIES

X = varargin{1};
if nargin < 2
    help iff
end
Y = varargin{2};

switch class(varargin{1})
    case {'lmi','constraint','sdpvar'}
        varargout{1} = setupMeta(lmi([]), mfilename,varargin{:});

    case 'char'
        varargout{1} = iff_internal(varargin{3:end});
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% MATLAB computation of pulse transformer model - DS Method
% File:  c:\M_files\shortcuts\xfrmrds2.m
% 9/19/02; 4/17/04; 2/15/07
%
tic;clc;clear;
K=1e3;pF=1e-12;mH=1e-3;uH=1e-6;ns=1e-9;ps=1e-12; % unit suffixes
%
% Components
%
R1=10;R2=1.5;R3=20*K;R4=1.5;R5=1*K;R6=0.5;R7=1;
C1=20*pF;C2=5*pF;C3=20*pF;L1=1*uH;L2=2*mH;L3=1*uH;
%
% Get A, B, D, & E arrays; this function called only once.
%
Nom=[R1 R2 R3 R4 R5 R6 R7 C1 C2 C3 L1 L2 L3];
[A,B,D,E,I]=tfrmr2(Nom);
%
% * * * * * * * * * * * * Frequency response * * * * * * * * * * * *
%
Ein=10; % Change Ein from 1V to 10V.
%
BF=2;ND=6;PD=50;NP=ND*PD+1;L=linspace(BF,BF+ND,NP);
%
% Since the output is vC3, we dont need the D and E arrays.
% The cv output below is [vC1 vC2 vC3 iL1 iL2 iL3]'
% (a column vector).  Hence we need vC3 or cv(3).
%
for i=1:NP
   F=10^L(i);s=2*pi*F*j;
   cx=(s*I-A)\B*Ein;
%   cy=D*cx+E*Ein; % cy not used
   Vo=abs(cx(3)); % vC3 = Vo
   Vf(i)=20*log10(Vo);
end
%
% * * * * * * * * * * * * * Transient response * * * * * * * * * * *
%
Tx=1/max(max(abs(A)));
disp('Shortest circuit time constant');Tx
%Per=input('Sweep time? (sec)');
% set Sweep time to 200ns = 200e-9 to match Spice run.
Per=200*ns;
kmax=1e5; % kmax increased due to fast time constant Tx
dt = 2*ps
N=6;
%dt=Per/kmax;N=6;
t1=linspace(0,Per,kmax);IV=zeros(N,kmax);
%
% input ramp parameters
%
p=Ein/(5*dt);b=6*dt;pw=5e4*dt;c=pw+6*dt;d=pw+11*dt;
Ea1=ramp1(p,t1(1),dt)-ramp1(p,t1(1),b)-ramp1(p,t1(1),c)+ramp1(p,t1(1),d);
% initialize k = 1
IV(:,1)=B*Ea1*dt;
%
% iterate for k = 2,3,...kmax
%
for k=2:kmax
   Eak=ramp1(p,t1(k),dt)-ramp1(p,t1(k),b)-ramp1(p,t1(k),c)+ramp1(p,t1(k),d);
   IV(:,k)=A*IV(:,k-1)*dt+B*Eak*dt+IV(:,k-1);
end
%
% Plot frequency response
%
subplot(2,1,1)
h=plot(L,Vf,'k');
set(h,'LineWidth',2);
grid on;
axis([BF BF+ND -40 30]);
XT=linspace(BF,BF+ND,7);
set(gca,'xtick',XT);
ylabel('dBV');title('AC Output Vc3');
xlabel('Log Freq(Hz)');
%
% Plot time response
%
subplot(2,1,2)
h=plot(t1/ns,IV(1,:),'k',t1/ns,IV(3,:),'r');
set(h,'LineWidth',2);
axis auto
grid on;ylabel('Volts');title('Transient response, Vc1 & Vc3');
xlabel('nsec');
legend('Vc1','Vc3');

figure(1) % display plot on screen.
%
disp(' ');disp('Execution time in seconds');
ET=toc


text:
function [pt_errs_phys, pts_moved_phys, TRE_phys, TREstd_phys] = move_points_tre_general(pts_mov_phys, ...
                        pts_fix_phys, Tr_phys, spc, lu_orig)
    if ~exist('lu_orig') || isempty(lu_orig)
        lu_orig = [0, 0, 0];
    end

    pts_moved_phys = 0 * pts_mov_phys;
    tsz = size(Tr_phys);
    tsz = tsz(1:3);
    for i = 1 : size(pts_mov_phys, 1)
        pos = round( (fl(pts_fix_phys(i, :))- lu_orig(:))./spc(:) );
        pos = min(pos, tsz(:));
        pos = max(pos, 1);
%                 pos
%         size(Tr_phys)
        pts_moved_phys(i, :) = (pts_fix_phys(i, :)' + fl(Tr_phys(pos(1), pos(2), pos(3), :)))';
    end

%     max(Tr_phys(:))
%     pts_moved_phys = bsxfun(@times, round(bsxfun(@rdivide, pts_moved_phys, spc)), spc);
%     pts_mov_phys = bsxfun(@times, round(bsxfun(@rdivide, pts_mov_phys, spc)), spc);

    pt_errs_phys = sqrt( sum((pts_moved_phys - pts_mov_phys).^2, 2) );

    TRE_phys = mean(pt_errs_phys);
    TREstd_phys = std(pt_errs_phys);
end

repo: pTVreg


text:
function J = pidtest(G,dt,parms)

s = tf('s');
K = parms(1) + parms(2)/s + parms(3)*s/(1+.001*s);
Loop = series(K,G);
ClosedLoop = feedback(Loop,1);
t = 0:dt:20;
[y,t] = step(ClosedLoop,t);

CTRLtf = K/(1+K*G);
u = lsim(K,1-y,t);

Q = 1;
R = .001;
J = dt*sum(Q*(1-y(:)).^2 + R*u(:).^2)

step(ClosedLoop,t)
h = findobj(gcf,'type','line');
set(h,'linewidth',2);
drawnow

author: dynamicslab


text:
function allen_ccf_npx_4shank(tv,av,st)
% allen_ccf_npx(tv,av,st)
% Andy Peters (peters.andrew.j@gmail.com)
% bunch of 4-shank modifications: Flora

% GUI for planning Neuropixels trajectories with the Allen CCF
% Part of repository: https://github.com/cortex-lab/allenCCF
% - directions for installing atlas in repository readme
% - some dependent functions from that repository
%
% (optional inputs - if CCF path written in Line 22, loaded automatically)
% tv, av, st = CCF template volume, annotated volume, structure tree

% Initialize gui_data structure
gui_data = struct;

% Allen CCF-bregma transform (estimated from eyeballing Paxinos->CCF)
% [AP,DV,ML]
bregma = [540,0,570];

% If not already loaded in, load in atlas
if nargin < 3
    allen_atlas_path = 'C:\Users\Andrew\OneDrive for Business\Documents\Atlases\AllenCCF';
    if isempty(allen_atlas_path)
        error('Enter path where Allen CCF is stored at Line 23');
    end
    tv = readNPY([allen_atlas_path filesep 'template_volume_10um.npy']); % grey-scale "background signal intensity"
    av = readNPY([allen_atlas_path filesep 'annotation_volume_10um_by_index.npy']); % the number at each pixel labels the area, see note below
    st = loadStructureTree([allen_atlas_path filesep 'structure_tree_safe_2017.csv']); % a table of what all the labels mean
end

% Load the colormap (located in the repository, find by associated fcn)
allenCCF_path = fileparts(which('allenCCFbregma'));
cmap_filename = [allenCCF_path filesep 'allen_ccf_colormap_2017.mat'];
load(cmap_filename);

% Set up the gui
probe_atlas_gui = figure('Toolbar','none','Menubar','none','color','w', ...
    'Name','Atlas-probe viewer','Units','normalized','Position',[0.2,0.2,0.7,0.7]);

% Set up the atlas axes
axes_atlas = subplot(1,5,1);
[~, brain_outline] = plotBrainGrid([],axes_atlas);
hold(axes_atlas,'on');
axis vis3d equal off manual
view([-30,25]);
caxis([0 300]);
[ap_max,dv_max,ml_max] = size(tv);
xlim([-10,ap_max+10])
ylim([-10,ml_max+10])
zlim([-10,dv_max+10])

% 4 shank parameters
shank_spacing=0:20:60;
shank_no=4;

for i=1:shank_no
% Set up the probe area axes                    % add 4 of these
    axes_probe_areas(i) = subplot(1,5,i+1);
    axes_probe_areas(i).ActivePositionProperty = 'position';
    set(axes_probe_areas(i),'FontSize',11);
    yyaxis(axes_probe_areas(i),'left');
    probe_areas_plot(i) = image(0);
    set(axes_probe_areas(i),'XTick','','YLim',[0,2880],'YColor','k','YDir','reverse');
    ylabel(axes_probe_areas(i),'');

    colormap(axes_probe_areas(i),cmap);
    caxis([1,size(cmap,1)])
    yyaxis(axes_probe_areas(i),'right');
    set(axes_probe_areas(i),'XTick','','YLim',[0,2880],'YColor','k','YDir','reverse');
    set(axes_probe_areas(i),'YTicklabels',[])
    title(axes_probe_areas(i),sprintf('shank %d',i));
end
yyaxis(axes_probe_areas(1),'left');
ylabel(axes_probe_areas(1),'Depth (\mum)');

% Position the axes              % change positioning so that there are 4
set(axes_atlas,'Position',[-0.15,-0.1,.8,1.2]);
set(axes_probe_areas(1),'Position',[.55,0.05,0.01,0.8]);
set(axes_probe_areas(2),'Position',[.67,0.05,0.01,0.8]);
set(axes_probe_areas(3),'Position',[.79,0.05,0.01,0.8]);
set(axes_probe_areas(4),'Position',[.91,0.05,0.01,0.8]);

% Set the current axes to the atlas (dirty, but some gca requirements)
axes(axes_atlas);

% Set up the probe reference/actual % again set up the reference so that 4
% shans are plotted
probe_ref_top = [bregma(1),bregma(3),0];
probe_ref_bottom = [bregma(1),bregma(3),size(tv,2)];
probe_ref_vector = [probe_ref_top',probe_ref_bottom'];
shank_spacing_vector=repmat(shank_spacing,2,1)';
allshanks_ref_vector=repmat(probe_ref_vector(2,:),shank_no,1)+shank_spacing_vector;

probe_ref_line = line(probe_ref_vector(1,:),allshanks_ref_vector,probe_ref_vector(3,:), ...
    'linewidth',1.5,'color','r','linestyle','--');

% probe_ref_line = line(probe_ref_vector(1,:),probe_ref_vector(2,:),probe_ref_vector(3,:), ...
%     'linewidth',1.5,'color','r','linestyle','--');

probe_length = 288.8; % IMEC phase 2.0 (in 10 ums)
probe_vector = [probe_ref_vector(:,1),diff(probe_ref_vector,[],2)./ ...
    norm(diff(probe_ref_vector,[],2))*probe_length + probe_ref_vector(:,1)];

allshanks_vector=repmat(probe_vector(2,:),shank_no,1)+shank_spacing_vector;
probe_line = line(probe_vector(1,:),allshanks_vector,probe_vector(3,:),'linewidth',3,'color','b','linestyle','-');
set(probe_line(1),'color','g'); % set shank 0 to be a different color

%probe_line = line(probe_vector(1,:),probe_vector(2,:),probe_vector(3,:),'linewidth',3,'color','b','linestyle','-');

% Set up the text to display coordinates
probe_coordinates_text = uicontrol('Style','text','String','', ...
    'Units','normalized','Position',[0,0.95,1,0.05], ...
    'BackgroundColor','w','HorizontalAlignment','left','FontSize',12);

% Store data
gui_data.tv = tv; % Intensity atlas
gui_data.av = av; % Annotated atlas
gui_data.st = st; % Labels table
gui_data.cmap = cmap; % Atlas colormap
gui_data.bregma = bregma; % Bregma for external referencing
gui_data.probe_length = probe_length; % Length of probe
gui_data.structure_plot_idx = []; % Plotted structures
gui_data.probe_angle = [0;0]; % Probe angles in ML/DV
gui_data.spherical_angle=[0;0];
gui_data.shank_no=shank_no;
gui_data.spin=0;
gui_data.shank_spacing=shank_spacing;

%Store handles
gui_data.handles.cortex_outline = brain_outline;
gui_data.handles.structure_patch = []; % Plotted structures
gui_data.handles.axes_atlas = axes_atlas; % Axes with 3D atlas
gui_data.handles.axes_probe_areas = axes_probe_areas; % Axes with probe areas
gui_data.handles.slice_plot = surface('EdgeColor','none'); % Slice on 3D atlas
gui_data.handles.slice_volume = 'tv'; % The volume shown in the slice
gui_data.handles.probe_ref_line = probe_ref_line; % Probe reference line on 3D atlas
gui_data.handles.probe_line = probe_line; % Probe reference line on 3D atlas
gui_data.handles.probe_areas_plot = probe_areas_plot; % Color-coded probe regions
gui_data.probe_coordinates_text = probe_coordinates_text; % Probe coordinates text

% Make 3D rotation the default state (toggle on/off with 'r')
h = rotate3d(axes_atlas);
h.Enable = 'on';
% Update the slice whenever a rotation is completed
h.ActionPostCallback = @update_slice;

% Set functions for key presses
hManager = uigetmodemanager(probe_atlas_gui);
[hManager.WindowListenerHandles.Enabled] = deal(false);
set(probe_atlas_gui,'KeyPressFcn',@key_press);
set(probe_atlas_gui,'KeyReleaseFcn',@key_release);

% Upload gui_data
guidata(probe_atlas_gui, gui_data);

% Display the first slice and update the probe position
update_slice(probe_atlas_gui);
update_probe_coordinates(probe_atlas_gui);

% Display controls
display_controls;

end

function key_press(probe_atlas_gui,eventdata)

% Get guidata
gui_data = guidata(probe_atlas_gui);
switch eventdata.Key

    case 'uparrow'
        if isempty(eventdata.Modifier)
            % Up: move probe anterior
            ap_offset = -10;

            for myshank=1:numel(gui_data.handles.probe_ref_line)
            set(gui_data.handles.probe_ref_line(myshank,1),'XData',get(gui_data.handles.probe_ref_line(myshank,1),'XData') + ap_offset);
            set(gui_data.handles.probe_line(myshank,1),'XData',get(gui_data.handles.probe_line(myshank,1),'XData') + ap_offset);
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-up: increase DV angle
            angle_change=-1;
            update_rotation(probe_atlas_gui,angle_change,'AP');
            % update gui value
%             new_angle = gui_data.probe_angle(2) + angle_change;
%             new_angle=check_for360(new_angle);
%             gui_data.probe_angle(2) = new_angle;

        elseif any(strcmp(eventdata.Modifier,'alt'))
            % Alt-up: raise probe
            probe_offset = -10;

            for myshank=1:numel(gui_data.handles.probe_line)
            old_probe_vector = cell2mat(get(gui_data.handles.probe_line(myshank,1),{'XData','YData','ZData'})');

            move_probe_vector = diff(old_probe_vector,[],2)./ ...
                norm(diff(old_probe_vector,[],2))*probe_offset;

            new_probe_vector = bsxfun(@plus,old_probe_vector,move_probe_vector);

            set(gui_data.handles.probe_line(myshank,1),'XData',new_probe_vector(1,:), ...
                'YData',new_probe_vector(2,:),'ZData',new_probe_vector(3,:));
            end
        end

    case 'downarrow'
        if isempty(eventdata.Modifier)
            % Down: move probe posterior
            ap_offset = 10;
            for myshank=1:numel(gui_data.handles.probe_ref_line)
            set(gui_data.handles.probe_ref_line(myshank,1),'XData',get(gui_data.handles.probe_ref_line(myshank,1),'XData') + ap_offset);
            set(gui_data.handles.probe_line(myshank,1),'XData',get(gui_data.handles.probe_line(myshank,1),'XData') + ap_offset);
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-down: decrease DV angle

            angle_change=1;
            update_rotation(probe_atlas_gui,angle_change,'AP');

%             new_angle = gui_data.probe_angle(2) + angle_change;
%             new_angle=check_for360(new_angle);
%             gui_data.probe_angle(2) = new_angle;

        elseif any(strcmp(eventdata.Modifier,'alt'))
            % Alt-down: lower probe
            probe_offset = 10;
            for myshank=1:numel(gui_data.handles.probe_line)  % this could be yoked together
            old_probe_vector = cell2mat(get(gui_data.handles.probe_line(myshank,1),{'XData','YData','ZData'})');

            move_probe_vector = diff(old_probe_vector,[],2)./ ...
                norm(diff(old_probe_vector,[],2))*probe_offset;

            new_probe_vector = bsxfun(@plus,old_probe_vector,move_probe_vector);

            set(gui_data.handles.probe_line(myshank,1),'XData',new_probe_vector(1,:), ...
                'YData',new_probe_vector(2,:),'ZData',new_probe_vector(3,:));
            end
        end

    case 'rightarrow'
        if isempty(eventdata.Modifier)
            % Right: move probe right
            ml_offset = 10;
            for myshank=1:numel(gui_data.handles.probe_ref_line)
            set(gui_data.handles.probe_ref_line(myshank,1),'YData',get(gui_data.handles.probe_ref_line(myshank,1),'YData') + ml_offset);
            set(gui_data.handles.probe_line(myshank,1),'YData',get(gui_data.handles.probe_line(myshank,1),'YData') + ml_offset);
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))

            angle_change=-1;
            update_rotation(probe_atlas_gui,angle_change,'ML');

%             new_angle = gui_data.probe_angle(1) + angle_change;
%             new_angle=check_for360(new_angle);
%
%             gui_data.probe_angle(1) = new_angle;
%

        elseif any(strcmp(eventdata.Modifier,'alt'))
            % alt-right : rotate clockwise
            angle_change=-1;
            update_rotation(probe_atlas_gui,angle_change,'spin');
            %update_spin(probe_atlas_gui,angle_change);
            new_angle=gui_data.spin+angle_change;
            new_angle=check_for360(new_angle);
            gui_data.spin=new_angle;

        end

    case 'leftarrow'
        if isempty(eventdata.Modifier)
            % Left: move probe left
            ml_offset = -10;
            for myshank=1:numel(gui_data.handles.probe_ref_line)
            set(gui_data.handles.probe_ref_line(myshank,1),'YData',get(gui_data.handles.probe_ref_line(myshank,1),'YData') + ml_offset);
            set(gui_data.handles.probe_line(myshank,1),'YData',get(gui_data.handles.probe_line(myshank,1),'YData') + ml_offset);
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))

            angle_change=1;
            update_rotation(probe_atlas_gui,angle_change,'ML');

%             new_angle = gui_data.probe_angle(1) + angle_change;
%             new_angle=check_for360(new_angle);
%             gui_data.probe_angle(1) = new_angle;

        elseif any(strcmp(eventdata.Modifier,'alt'))
            % alt-left : rotate anticlockwise
            angle_change=1;
            update_rotation(probe_atlas_gui,angle_change,'spin');
            %update_spin(probe_atlas_gui,angle_change');
            new_angle=gui_data.spin+angle_change;
            new_angle=check_for360(new_angle);
            gui_data.spin=new_angle;
        end

    case 'c'
        % Bring up controls again
        display_controls;

    case 'b'
        % Toggle brain outline visibility
        current_visibility = gui_data.handles.cortex_outline.Visible;
        switch current_visibility; case 'on'; new_visibility = 'off'; case 'off'; new_visibility = 'on'; end;
        set(gui_data.handles.cortex_outline,'Visible',new_visibility);

    case 'a'
        % Toggle plotted structure visibility
        if ~isempty(gui_data.structure_plot_idx)
            current_visibility = get(gui_data.handles.structure_patch(1),'Visible');
            switch current_visibility; case 'on'; new_visibility = 'off'; case 'off'; new_visibility = 'on'; end;
            set(gui_data.handles.structure_patch,'Visible',new_visibility);
        end

    case 's'
        % Toggle slice volume/visibility
        slice_volumes = {'tv','av','none'};
        new_slice_volume = slice_volumes{circshift( ...
            strcmp(gui_data.handles.slice_volume,slice_volumes),[0,1])};

        if strcmp(new_slice_volume,'none')
            set(gui_data.handles.slice_plot,'Visible','off');
        else
            set(gui_data.handles.slice_plot,'Visible','on');
        end

        gui_data.handles.slice_volume = new_slice_volume;
        guidata(probe_atlas_gui, gui_data);

        update_slice(probe_atlas_gui);

    case 'p'
        % Toggle probe visibility
        current_visibility = gui_data.handles.probe_ref_line.Visible;
        switch current_visibility; case 'on'; new_visibility = 'off'; case 'off'; new_visibility = 'on'; end;
        set(gui_data.handles.probe_ref_line,'Visible',new_visibility);
        set(gui_data.handles.probe_line,'Visible',new_visibility);

    case 'r'
        % Toggle 3D rotation
        h = rotate3d(gui_data.handles.axes_atlas);
        switch h.Enable
            case 'off'
                h.Enable = 'on';
                % Update the slice whenever a rotation is completed
                h.ActionPostCallback = @update_slice;
                %(need to restore key-press functionality with rotation)
                hManager = uigetmodemanager(probe_atlas_gui);
                [hManager.WindowListenerHandles.Enabled] = deal(false);
                set(probe_atlas_gui,'KeyPressFcn',@key_press);
            case 'on'
                h.Enable = 'off';
        end

    case 'm'
        % Set probe angle
        set_probe_position(probe_atlas_gui);
        % Get updated guidata
        gui_data = guidata(probe_atlas_gui);

    case {'equal','add'}
        % Add structure(s) to display
        slice_spacing = 10;

        % Prompt for which structures to show (only structures which are
        % labelled in the slice-spacing downsampled annotated volume)

        if ~any(strcmp(eventdata.Modifier,'shift'))
            % (no shift: list in native CCF order)
            parsed_structures = unique(reshape(gui_data.av(1:slice_spacing:end, ...
                1:slice_spacing:end,1:slice_spacing:end),[],1));
            plot_structures_parsed = listdlg('PromptString','Select a structure to plot:', ...
                'ListString',gui_data.st.safe_name(parsed_structures),'ListSize',[520,500]);
            plot_structures = parsed_structures(plot_structures_parsed);

            if ~isempty(plot_structures)
                for curr_plot_structure = reshape(plot_structures,1,[])
                    % If this label isn't used, don't plot
                    if ~any(reshape(gui_data.av( ...
                            1:slice_spacing:end,1:slice_spacing:end,1:slice_spacing:end),[],1) == curr_plot_structure)
                        disp(['"' gui_data.st.safe_name{curr_plot_structure} '" is not parsed in the atlas'])
                        continue
                    end

                    gui_data.structure_plot_idx(end+1) = curr_plot_structure;

                    plot_structure_color = hex2dec(reshape(gui_data.st.color_hex_triplet{curr_plot_structure},2,[])')./255;
                    structure_3d = isosurface(permute(gui_data.av(1:slice_spacing:end, ...
                        1:slice_spacing:end,1:slice_spacing:end) == curr_plot_structure,[3,1,2]),0);

                    structure_alpha = 0.2;
                    gui_data.handles.structure_patch(end+1) = patch('Vertices',structure_3d.vertices*slice_spacing, ...
                        'Faces',structure_3d.faces, ...
                        'FaceColor',plot_structure_color,'EdgeColor','none','FaceAlpha',structure_alpha);
                end
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % (shift: use hierarchy search)
            plot_structures = hierarchicalSelect(gui_data.st);

            if ~isempty(plot_structures) % will be empty if dialog was cancelled
                % get all children of this one
                thisID = gui_data.st.id(plot_structures);
                idStr = sprintf('/%d/', thisID);
                theseCh = find(cellfun(@(x)contains(x,idStr), gui_data.st.structure_id_path));

                % plot the structure
                slice_spacing = 5;
                plot_structure_color = hex2dec(reshape(gui_data.st.color_hex_triplet{plot_structures},3,[]))./255;
                structure_3d = isosurface(permute(ismember(gui_data.av(1:slice_spacing:end, ...
                    1:slice_spacing:end,1:slice_spacing:end),theseCh),[3,1,2]),0);

                structure_alpha = 0.2;
                gui_data.structure_plot_idx(end+1) = plot_structures;
                gui_data.handles.structure_patch(end+1) = patch('Vertices',structure_3d.vertices*slice_spacing, ...
                    'Faces',structure_3d.faces, ...
                    'FaceColor',plot_structure_color,'EdgeColor','none','FaceAlpha',structure_alpha);

            end

        end

    case {'hyphen','subtract'}
        % Remove structure(s) already plotted
        if ~isempty(gui_data.structure_plot_idx)
            remove_structures = listdlg('PromptString','Select a structure to remove:', ...
                'ListString',gui_data.st.safe_name(gui_data.structure_plot_idx));
            delete(gui_data.handles.structure_patch(remove_structures))
            gui_data.structure_plot_idx(remove_structures) = [];
            gui_data.handles.structure_patch(remove_structures) = [];
        end

    case 'x'
        % Export the probe coordinates in Allen CCF to the workspace
        probe_vector = cell2mat(get(gui_data.handles.probe_line,{'XData','YData','ZData'})');
        probe_vector_ccf = round(probe_vector([1,3,2],:))';
        assignin('base','probe_vector_ccf',probe_vector_ccf)
        disp('Copied probe vector coordinates to workspace');

    case 'h'
        % Load probe histology points, plot line of best fit
        probe_histology(probe_atlas_gui);
end

% Upload gui_data
guidata(probe_atlas_gui, gui_data);

end

function key_release(probe_atlas_gui,eventdata)

% Get guidata
gui_data = guidata(probe_atlas_gui);

switch eventdata.Key
    case {'rightarrow','leftarrow','uparrow','downarrow'}
        % Update the probe info/slice on arrow release
        update_probe_coordinates(probe_atlas_gui);
        update_slice(probe_atlas_gui);
end

% Upload gui_data
guidata(probe_atlas_gui, gui_data);

end

function update_slice(probe_atlas_gui,varargin)

% Get guidata
gui_data = guidata(probe_atlas_gui);

% Only update the slice if it's visible
if strcmp(gui_data.handles.slice_plot(1).Visible,'on')

    % Get current position of camera
    curr_campos = campos;

    % Get probe vector
    probe_ref_top = [gui_data.handles.probe_ref_line(1,1).XData(1), ...
        gui_data.handles.probe_ref_line(1,1).YData(1),gui_data.handles.probe_ref_line(1,1).ZData(1)];
    probe_ref_bottom = [gui_data.handles.probe_ref_line(1).XData(2), ...
        gui_data.handles.probe_ref_line(1,1).YData(2),gui_data.handles.probe_ref_line(1,1).ZData(2)];
    probe_vector = probe_ref_top - probe_ref_bottom;

    % Get probe-camera vector
    probe_camera_vector = probe_ref_top - curr_campos;

    % Get the vector to plot the plane in (along with probe vector)
    plot_vector = cross(probe_camera_vector,probe_vector);

    % Get the normal vector of the plane
    normal_vector = cross(plot_vector,probe_vector);

    % Get the plane offset through the probe
    plane_offset = -(normal_vector*probe_ref_top');

    % Define a plane of points to index
    % (the plane grid is defined based on the which cardinal plan is most
    % orthogonal to the plotted plane. this is janky but it works)
    slice_px_space = 3;
    %[~,cam_plane] = max(abs((campos - camtarget)./norm(campos - camtarget)));

    [~,cam_plane] = max(abs(normal_vector./norm(normal_vector)));

    switch cam_plane

        case 1
            [plane_y,plane_z] = meshgrid(1:slice_px_space:size(gui_data.tv,3),1:slice_px_space:size(gui_data.tv,2));
            plane_x = ...
                (normal_vector(2)*plane_y+normal_vector(3)*plane_z + plane_offset)/ ...
                -normal_vector(1);

        case 2
            [plane_x,plane_z] = meshgrid(1:slice_px_space:size(gui_data.tv,1),1:slice_px_space:size(gui_data.tv,2));
            plane_y = ...
                (normal_vector(1)*plane_x+normal_vector(3)*plane_z + plane_offset)/ ...
                -normal_vector(2);

        case 3
            [plane_x,plane_y] = meshgrid(1:slice_px_space:size(gui_data.tv,1),1:slice_px_space:size(gui_data.tv,3));
            plane_z = ...
                (normal_vector(1)*plane_x+normal_vector(2)*plane_y + plane_offset)/ ...
                -normal_vector(3);

    end

    % Get the coordiates on the plane
    x_idx = round(plane_x);
    y_idx = round(plane_y);
    z_idx = round(plane_z);

    % Find plane coordinates in bounds with the volume
    use_xd = x_idx > 0 & x_idx < size(gui_data.tv,1);
    use_yd = y_idx > 0 & y_idx < size(gui_data.tv,3);
    use_zd = z_idx > 0 & z_idx < size(gui_data.tv,2);
    use_idx = use_xd & use_yd & use_zd;

    curr_slice_idx = sub2ind(size(gui_data.tv),x_idx(use_idx),z_idx(use_idx),y_idx(use_idx));

    % Find plane coordinates that contain brain
    curr_slice_isbrain = false(size(use_idx));
    curr_slice_isbrain(use_idx) = gui_data.av(curr_slice_idx) > 1;

    % Index coordinates in bounds + with brain
    grab_pix_idx = sub2ind(size(gui_data.tv),x_idx(curr_slice_isbrain),z_idx(curr_slice_isbrain),y_idx(curr_slice_isbrain));

    % Grab pixels from (selected) volume
    curr_slice = nan(size(use_idx));
    switch gui_data.handles.slice_volume
        case 'tv'
            curr_slice(curr_slice_isbrain) = gui_data.tv(grab_pix_idx);
            colormap(gui_data.handles.axes_atlas,'gray');
            caxis([0,255]);
        case 'av'
            curr_slice(curr_slice_isbrain) = gui_data.av(grab_pix_idx);
            colormap(gui_data.handles.axes_atlas,gui_data.cmap);
            caxis([1,size(gui_data.cmap,1)]);
    end

    % Update the slice display
    set(gui_data.handles.slice_plot,'XData',plane_x,'YData',plane_y,'ZData',plane_z,'CData',curr_slice);

    % Upload gui_data
    guidata(probe_atlas_gui, gui_data);

end

end

function set_probe_position(probe_atlas_gui,varargin)

% Get guidata
gui_data = guidata(probe_atlas_gui);

% Prompt for angles
prompt_text = { ...
    'AP position (mm from bregma)', ...
    'ML position (mm from bregma)', ...
    'ML angle (relative to midline)', ...
    'DV angle (relative to horizontal)',...
    'probe spin'};

new_probe_position = cellfun(@str2num,inputdlg(prompt_text,'Set probe position',1));

% Convert probe position: mm->CCF and degrees->radians
probe_ccf_coordinates = round(gui_data.bregma([1,3])' - new_probe_position(1:2)*100);
%probe_angle_rad = (new_probe_position(3:4)/360)*2*pi;

% Update the probe and trajectory reference
[ap_max,dv_max,ml_max] = size(gui_data.tv);

max_ref_length = sqrt(sum(([ap_max,dv_max,ml_max].^2)));

probe_angle_rad = (new_probe_position(3:4)/360)*2*pi;
[x,y,z] = sph2cart(probe_angle_rad(1),probe_angle_rad(2),max_ref_length);

% Get top of probe reference with user brain intersection point
% (get DV location of brain surface at point)
probe_brain_dv = find(gui_data.av(probe_ccf_coordinates(1),:, ...
    probe_ccf_coordinates(2)) > 1,1);
% (back up to 0 DV in CCF space)
probe_ref_top_ap = interp1(probe_brain_dv+[0,z],probe_ccf_coordinates(1)+[0,x],0,'linear','extrap');
probe_ref_top_ml = interp1(probe_brain_dv+[0,z],probe_ccf_coordinates(2)+[0,y],0,'linear','extrap');

% Set new probe position
probe_ref_top = [probe_ref_top_ap,probe_ref_top_ml,0];
probe_ref_bottom = probe_ref_top + [x,y,z];
probe_ref_vector = [probe_ref_top;probe_ref_bottom]';

%[~,,~] = sph2cart(probe_angle_rad(2),probe_angle_rad(1),max_ref_length);

% probe_ref_top = [probe_ccf_coordinates(1),probe_ccf_coordinates(2),0];
% probe_ref_bottom = [probe_ccf_coordinates(1),probe_ccf_coordinates(2),max_ref_length];
% probe_ref_vector = [probe_ref_top;probe_ref_bottom]';

probe_vector = [probe_ref_vector(:,1),diff(probe_ref_vector,[],2)./ ...
    norm(diff(probe_ref_vector,[],2))*gui_data.probe_length + probe_ref_vector(:,1)];

shank_spacing_vector=reshape([gui_data.shank_spacing;gui_data.shank_spacing],[],1)';

probe_ref_vector=repmat(probe_ref_vector,1,gui_data.shank_no);
probe_ref_vector(2,:)=probe_ref_vector(2,:)+shank_spacing_vector;

probe_vector=repmat(probe_vector,1,gui_data.shank_no);
probe_vector(2,:)=probe_vector(2,:)+shank_spacing_vector;

gui_data.probe_angle(1) =new_probe_position(3);
gui_data.probe_angle(2) =new_probe_position(4);
gui_data.spin=new_probe_position(5);

probe_top=probe_vector(:,1);
probe_bottom=probe_vector(:,2);
u=(probe_bottom-probe_top)/norm(probe_bottom-probe_top);

rotation_matrix=get_rotation_matrix(gui_data.spin,'spin',u);
[probe_ref_vector,probe_vector]=please_rotate(probe_ref_vector,probe_vector,rotation_matrix);

for myshank=1:gui_data.shank_no
    st_ix=2*myshank-1;
    set(gui_data.handles.probe_ref_line(myshank),'XData',probe_ref_vector(1,st_ix:(st_ix+1)), ...
            'YData',probe_ref_vector(2,st_ix:(st_ix+1)), ...
            'ZData',probe_ref_vector(3,st_ix:(st_ix+1)));

    set(gui_data.handles.probe_line(myshank),'XData',probe_vector(1,st_ix:(st_ix+1)), ...
        'YData',probe_vector(2,st_ix:(st_ix+1)), ...
        'ZData',probe_vector(3,st_ix:(st_ix+1)));
end

% Upload gui_data
guidata(probe_atlas_gui, gui_data);

% Update the slice and probe coordinates
update_slice(probe_atlas_gui);
update_probe_coordinates(probe_atlas_gui);

end

function update_probe_coordinates(probe_atlas_gui,varargin)

% Get guidata
gui_data = guidata(probe_atlas_gui);

for myshank=1:numel(gui_data.handles.probe_ref_line)
    % Get the positions of the probe and trajectory reference
    probe_ref_vector = cell2mat(get(gui_data.handles.probe_ref_line(myshank,1),{'XData','YData','ZData'})');
    probe_vector = cell2mat(get(gui_data.handles.probe_line(myshank,1),{'XData','YData','ZData'})');

    % try reading the spherical coordinate here
    [theta,rho]=get_spherical_coords(probe_vector);

    %%
    trajectory_n_coords = max(abs(diff(probe_ref_vector,[],2)));
    [trajectory_xcoords,trajectory_ycoords,trajectory_zcoords] = deal( ...
        linspace(probe_ref_vector(1,1),probe_ref_vector(1,2),trajectory_n_coords), ...
        linspace(probe_ref_vector(2,1),probe_ref_vector(2,2),trajectory_n_coords), ...
        linspace(probe_ref_vector(3,1),probe_ref_vector(3,2),trajectory_n_coords));

    probe_n_coords = sqrt(sum(diff(probe_vector,[],2).^2));
    [probe_xcoords,probe_ycoords,probe_zcoords] = deal( ...
        linspace(probe_vector(1,1),probe_vector(1,2),probe_n_coords), ...
        linspace(probe_vector(2,1),probe_vector(2,2),probe_n_coords), ...
        linspace(probe_vector(3,1),probe_vector(3,2),probe_n_coords));

    % Get brain labels across the probe and trajectory, and intersection with brain
    pixel_space = 5;
    trajectory_areas = interp3(single(gui_data.av(1:pixel_space:end,1:pixel_space:end,1:pixel_space:end)), ...
        round(trajectory_zcoords/pixel_space),round(trajectory_xcoords/pixel_space),round(trajectory_ycoords/pixel_space),'nearest');
    trajectory_brain_idx = find(trajectory_areas > 1,1);
    trajectory_brain_intersect = ...
        [trajectory_xcoords(trajectory_brain_idx),trajectory_ycoords(trajectory_brain_idx),trajectory_zcoords(trajectory_brain_idx)]';

    probe_areas = interp3(single(gui_data.av(1:pixel_space:end,1:pixel_space:end,1:pixel_space:end)), ...
        round(probe_zcoords/pixel_space),round(probe_xcoords/pixel_space),round(probe_ycoords/pixel_space),'nearest')';
    probe_area_boundaries = intersect(unique([find(~isnan(probe_areas),1,'first'); ...
        find(diff(probe_areas) ~= 0);find(~isnan(probe_areas),1,'last')]),find(~isnan(probe_areas)));
    probe_area_centers = probe_area_boundaries(1:end-1) + diff(probe_area_boundaries)/2;
    probe_area_labels = gui_data.st.safe_name(probe_areas(round(probe_area_centers)));

    % Get position of brain intersect relative to bregma
    probe_bregma_coordinate = round((gui_data.bregma([1,3])' - trajectory_brain_intersect(1:2))*10);

    % Get the depth of the bottom of the probe (sign: hack by z offset)
    probe_depth = round(sqrt(sum((trajectory_brain_intersect - probe_vector(:,2)).^2))*10)* ...
        sign(probe_vector(3,2)-trajectory_brain_intersect(3));

    % Update the text
    if myshank==1
        gui_data.probe_angle(1)=rho;
        gui_data.probe_angle(2)=theta;

        probe_text = ['Shank 1 insertion: ' ...
                num2str(probe_bregma_coordinate(1)) ' AP, ', ...
                num2str(-probe_bregma_coordinate(2)) ' ML, ', ...
                num2str(probe_depth) ' Probe-axis, ', ...
                num2str(round(theta)) char(176) ' from midline, ' ...
                num2str(round(rho)) char(176) ' from horizontal, '  ...
                num2str(gui_data.spin) char(176) ' spin '];

       %         num2str(gui_data.probe_angle(2)) char(176) ' from midline, ' ...
        %        num2str(gui_data.probe_angle(1)) char(176) ' from horizontal, ' ...

        set(gui_data.probe_coordinates_text,'String',probe_text);
    end

    % Update the probe areas
    yyaxis(gui_data.handles.axes_probe_areas(myshank),'right');
    set(gui_data.handles.probe_areas_plot(myshank),'YData',[1:length(probe_areas)]*10,'CData',probe_areas);
    set(gui_data.handles.axes_probe_areas(myshank),'YTick',probe_area_centers*10,'YTickLabels',probe_area_labels','YTickLabelRotation',45);

end
% Upload gui_data
guidata(probe_atlas_gui, gui_data);

end

function probe_histology(probe_atlas_gui)
% Load histology points
% UNDER CONSTRUCTION
% (used to use SHARP-Track, now using mine)

% Get guidata
gui_data = guidata(probe_atlas_gui);

[probe_file,probe_path] = uigetfile('*.mat','Choose probe coordinate file');
load([probe_path,probe_file]);

if exist('pointList','var')
    histology_points = pointList.pointList{1};
elseif exist('probe_ccf','var')
    histology_points = probe_ccf(1).points; % only use first probe
end

r0 = mean(histology_points,1);
xyz = bsxfun(@minus,histology_points,r0);
[~,~,V] = svd(xyz,0);
histology_probe_direction = V(:,1);

probe_eval_points = [-1000,1000];
probe_line_endpoints = bsxfun(@plus,bsxfun(@times,probe_eval_points',histology_probe_direction'),r0);

% Philip's GUI: not saved in native CCF order?
% plot3(histology_points(:,3),histology_points(:,1),histology_points(:,2),'.b','MarkerSize',20);
% line(P(:,3),P(:,1),P(:,2),'color','k','linewidth',2)

% % Mine: saved in native CCF order [AP,DV,ML]
plot3(histology_points(:,1),histology_points(:,3),histology_points(:,2),'.b','MarkerSize',20);
line(probe_line_endpoints(:,1),probe_line_endpoints(:,3),probe_line_endpoints(:,2),'color','k','linewidth',2)

% Place the probe on the histology best-fit axis
[ap_max,dv_max,ml_max] = size(gui_data.tv);

probe_ref_top = probe_line_endpoints(1,[1,3,2]);
probe_ref_bottom = probe_line_endpoints(2,[1,3,2]);
probe_ref_vector = [probe_ref_top;probe_ref_bottom]';

set(gui_data.handles.probe_ref_line,'XData',probe_ref_vector(1,:), ...
    'YData',probe_ref_vector(2,:), ...
    'ZData',probe_ref_vector(3,:));

probe_vector = [probe_ref_vector(:,1),diff(probe_ref_vector,[],2)./ ...
    norm(diff(probe_ref_vector,[],2))*gui_data.probe_length + probe_ref_vector(:,1)];
set(gui_data.handles.probe_line,'XData',probe_vector(1,:), ...
    'YData',probe_vector(2,:),'ZData',probe_vector(3,:));

% Upload gui_data
[theta,phi] = cart2sph(diff(probe_ref_vector(1,:)),diff(probe_ref_vector(2,:)),diff(probe_ref_vector(3,:)));
gui_data.probe_angle = ([theta,phi]/(2*pi))*360;

guidata(probe_atlas_gui, gui_data);

% Update the slice and probe coordinates
update_slice(probe_atlas_gui);
update_probe_coordinates(probe_atlas_gui);

end

function display_controls

% Print controls
CreateStruct.Interpreter = 'tex';
CreateStruct.WindowStyle = 'non-modal';
msgbox( ...
    {'\fontsize{12}' ...
    '\bf Probe: \rm' ...
    'Arrow keys : translate probe' ...
    'Alt/Option up/down : raise/lower probe' ...
    'Alt/Option left/right : spin probe' ...
    'Shift arrow keys : change probe angle' ...
    'm : set probe location manually', ...
    '\bf 3D brain areas: \rm' ...
    ' =/+ : add (list selector)' ...
    ' Shift =/+ : add (hierarchy selector)' ...
    ' - : remove', ...
    '\bf Visibility: \rm' ...
    's : atlas slice (toggle tv/av/off)' ...
    'b : brain outline' ...
    'p : probe' ...
    'a : 3D brain areas' ...
    '\bf Other: \rm' ...
    'r : toggle clickable rotation' ...
    'x : export probe coordinates to workspace' ...
    'h : load and plot histology-defined trajectory', ...
    'c : bring up controls box'}, ...
    'Controls',CreateStruct);

end

function update_rotation(probe_atlas_gui,angle_change,mydir)
% this is the function when we update the rotation of the probe either in
% ML or AP
    gui_data = guidata(probe_atlas_gui);
    probe_ref_vector = cell2mat(get(gui_data.handles.probe_ref_line,{'XData','YData','ZData'})');
    probe_vector = cell2mat(get(gui_data.handles.probe_line,{'XData','YData','ZData'})');

    probe_top=probe_vector(:,1);
    probe_bottom=probe_vector(:,2);
    u=(probe_bottom-probe_top)/norm(probe_bottom-probe_top);

    rotation_matrix=get_rotation_matrix(angle_change,mydir,u);
    [rotated_ref_bregma,rotated_probe_bregma]=please_rotate(probe_ref_vector,probe_vector,rotation_matrix);

    for myshank=1:gui_data.shank_no
        st_ix=2*myshank-1;
        set(gui_data.handles.probe_ref_line(myshank),'XData',rotated_ref_bregma(1,st_ix:(st_ix+1)), ...
                'YData',rotated_ref_bregma(2,st_ix:(st_ix+1)), ...
                'ZData',rotated_ref_bregma(3,st_ix:(st_ix+1)));

        set(gui_data.handles.probe_line(myshank),'XData',rotated_probe_bregma(1,st_ix:(st_ix+1)), ...
            'YData',rotated_probe_bregma(2,st_ix:(st_ix+1)), ...
            'ZData',rotated_probe_bregma(3,st_ix:(st_ix+1)));
    end
end

function [rotation_matrix]=get_rotation_matrix(angle,direction,varargin)
angle_rad=deg2rad(angle);
    switch direction
        case 'ML'
            u=[1,0,0];
        case 'AP'
            u=[0,1,0];
        case 'spin' % that is any arbitrary unit vector that one requests
            u=varargin{1};
    end
    dotU=u*u';
    cross_M=[0,-u(3),u(2);u(3),0,-u(1);-u(2),u(1),0];
    rotation_matrix=cos(angle_rad)*eye(3)+sin(angle_rad)*cross_M+(1-cos(angle_rad))*dotU;
    rotation_matrix=round(rotation_matrix,3);
end

function [rotated_ref_bregma,rotated_probe_bregma]=please_rotate(probe_ref_vector,probe_vector,rotation_matrix)

    zeropoint=probe_vector(:,1);
    zerocentered_ref=probe_ref_vector-zeropoint;

    zerocentered_probe=probe_vector-zeropoint;

    rotated_ref=rotation_matrix*zerocentered_ref;
    rotated_ref_bregma=rotated_ref(1:3,:)+zeropoint;

    rotated_probe=rotation_matrix*zerocentered_probe;
    rotated_probe_bregma=rotated_probe(1:3,:)+zeropoint;
end

function [angle]=check_for360(angle)
    if angle>360
        angle=angle-360;
    elseif angle<-360
        angle=angle+360;
    end
end

function [theta,rho]=get_spherical_coords(probe_ref_vector)

[probe_azimuth,probe_elevation] = cart2sph( ...
    diff(fliplr(probe_ref_vector(1,1:2))), ...
    diff(fliplr(probe_ref_vector(2,1:2))), ...
    diff(fliplr(probe_ref_vector(3,1:2))));

a = [probe_azimuth,probe_elevation]*(360/(2*pi));
theta=a(1);
rho=a(2);

end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function [x, y] = findLinePoints(p1,p2);
%
%   [x, y] = findLinePoints(p1,p2)
%
% Author:  I think it was me.
% Purpose:
%      Find the x,y values that fall along a line between p1 and p2.
%
% 2002.07.24 RFD & FWC- fixed bug in vertical and horizontal special cases.
%            (Was returning column vectors instead of row.)

x1 = p1(1); y1 = p1(2);
x2 = p2(1); y2 = p2(2);

if y2 == y1
  if x1 == x2
    error;
    return;
  end
  x = [x1:x2];  y = y1*ones(1,length(x));
elseif x1 == x2
  if y1 == y2
    error;
    return;
  end
  y = [y1:y2];  x = x1*ones(1,length(y));
else
  slope = (y2-y1)/(x2-x1);
  b = y1 - slope*x1;
  if abs(y2 - y1) > abs(x2 - x1)
    if y1 < y2,  y = y1:y2;
    else,        y = y2:y1;
    end
    x = round( (y - b) / slope);
  else
    if x1 < x2, x = x1:x2;
    else        x = x2:x1;
    end
    y = round(slope*x + b);
  end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%
% Compute features for a set of video files from datasets
%
close all;
clear;
warning('off','all');
addpath(genpath('./lib/iqm/FRIQUEE_feat_sel_release'));

%%
% parameters
algo_name = 'FRIQUEE_feat_sel';
data_name = 'KONVID_1K';

%% *You need to customize here*
if strcmp(data_name, 'TEST_VIDEOS')
    data_path = 'videos'; % dataset video path
elseif strcmp(data_name, 'KONVID_1K')
    data_path = '/media/ztu/Seagate-ztu-ugc/KONVID_1K/KoNViD_1k_videos';
elseif strcmp(data_name, 'LIVE_VQC')
    data_path = '/media/ztu/Seagate-ztu-ugc/LIVE_VQC/VideoDatabase';
elseif strcmp(data_name, 'YOUTUBE_UGC')
    data_path = '/media/ztu/Seagate-ztu-ugc/YT_UGC/original_videos';
elseif strcmp(data_name, 'LIVE_VQA')
    data_path = '/media/ztu/Seagate-ztu/LIVE_VQA/videos';
end

video_tmp = 'video_tmp';
if ~exist(video_tmp, 'dir'), mkdir(video_tmp); end
feat_path = '../features';
filelist_csv = fullfile(feat_path, [data_name,'_metadata.csv']);
filelist = readtable(filelist_csv);
num_videos = size(filelist,1);
out_path = './feat_sel_mats';
if ~exist(out_path, 'dir'), mkdir(out_path); end
out_feat_name = fullfile(out_path, [data_name,'_',algo_name,'_feats.mat']);
out_feat_frames_name = fullfile(out_path, [data_name,'_',algo_name,'_frames_feats.mat']);
feats_mat = [];
feats_mat_frames = cell( num_videos, 1 );
%===================================================

tic
for i = 1:num_videos

    % get video full path and decoded video name
	if strcmp(data_name, 'TEST_VIDEOS')
	    video_name = fullfile(data_path,  filelist.video_name{i});
	    yuv_name = fullfile(video_tmp, [filelist.video_name{i}, '.yuv']);
	elseif strcmp(data_name, 'KONVID_1K')
	    video_name = fullfile(data_path,  [num2str(filelist.flickr_id(i)),'.mp4']);
	    yuv_name = fullfile(video_tmp, [num2str(filelist.flickr_id(i)), '.yuv']);
	elseif strcmp(data_name, 'LIVE_VQC')
	    video_name = fullfile(data_path, filelist.File{i});
	    yuv_name = fullfile(video_tmp, [filelist.File{i}, '.yuv']);
	elseif strcmp(data_name, 'YOUTUBE_UGC')
	    video_name = fullfile(data_path, filelist.category{i},...
		[num2str(filelist.resolution(i)),'P'],[filelist.vid{i},'.mkv']);
	    yuv_name = fullfile(video_tmp, [filelist.vid{i}, '.yuv']);
	elseif strcmp(data_name, 'LIVE_VQA')
	    strs = strsplit(filelist.filename{i}, '_');
	    video_name = fullfile(data_path, [strs{1}(1:2), '_Folder'], filelist.filename{i});
	    yuv_name = video_name;
	end
	fprintf('\n---\nComputing features for %d-th sequence: %s\n', i, video_name);
    if ~strcmp(video_name, yuv_name)
    cmd = ['ffmpeg -loglevel error -y -i ', video_name, ' -pix_fmt yuv420p -vsync 0 ', yuv_name];
    system(cmd);
    end

    % get video meta data
    width = filelist.width(i);
    height = filelist.height(i);
    framerate = round(filelist.framerate(i));
    nb_frames = filelist.nb_frames(i);

    % read YUV frame (credit: Dae Yeol Lee)
    fp_input = fopen(yuv_name, 'r');
    uv_width = width/2;
    uv_height = height/2;
    feats_frames = [];

    tic
    for fr = 1:2:nb_frames-5
%         fr % frame No. printer

        try
            %% Start a file pointer
            fseek(fp_input,(fr-1)*1.5*width*height, 'bof'); % Frame read for 8 bit
            %% Y component
            %1) read y stream
            y_stream = fread(fp_input, width * height, 'uchar'); % for 8 bit
            % 2) reshape into a plane
            y_plane= reshape(y_stream, width, height).';
            %% U component
            %1) read u stream
            u_stream = fread(fp_input, uv_width * uv_height, 'uchar'); % for 8 bit
            % 2) reshape into a plane
            u_plane= reshape(u_stream, uv_width, uv_height).';
            %% V component
            %1) read v stream
            v_stream = fread(fp_input, uv_width * uv_height, 'uchar'); % for 8 bit
            % 2) reshape into a plane
            v_plane= reshape(v_stream, uv_width, uv_height).';
            % yuv2rgb
            u_plane = imresize(u_plane, size(y_plane), 'bicubic');
            v_plane = imresize(v_plane, size(y_plane), 'bicubic');
            rgb_plane = ycbcr2rgb(uint8(cat(3, y_plane, u_plane, v_plane)));

            % extract features
%             tic
            feat_tmp = extractFRIQUEEFeatures(rgb_plane);
%             toc
            feats_frames(end+1,:) = feat_tmp.friqueeALL;

        catch
            continue
        end

    end
    fclose(fp_input);
    delete(yuv_name)
    feats_mat_frames{i} = feats_frames;

        % compute brisque mean and consistency within each 1-sec chunk!!
    cons_feats = [];
    mean_feats = [];
    n_temp_vecs = length(feats_frames(:,1));
    half_blk_len = floor(framerate/2);

    fprintf('Pooling mean and consistency features\n');

    for j = 1:half_blk_len:n_temp_vecs - half_blk_len

        j_start = j;
        j_end = j + half_blk_len;

        % compute consistency features
        cons_feats = [cons_feats; std(feats_frames(j_start:j_end, :))];
        mean_feats = [mean_feats; mean(feats_frames(j_start:j_end, :))];
    end

    feats = [mean(mean_feats) ...
             mean(cons_feats) ];

    feats_mat(i,:) = feats;
    toc
end
toc
save(out_feat_name, 'feats_mat');
save(out_feat_frames_name, 'feats_mat_frames');

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

% Find Indices of Positive and Negative Examples

pos = find(y==1);
neg = find(y == 0);
% Plot Examples

plot(X(pos, 1), X(pos, 2), 'k+','LineWidth', 2,'MarkerSize', 7);
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

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

hold off;

end

author: xjwhhh

repo: AndrewNgMachineLearning


text:
% File Name: resampling.m
% Author: Heung-Il Suk
% Cite: H.-I. Suk and S.-W. Lee, "A Novel Bayesian Framework for Discriminative
% Feature Extraction in Brain-Computer Interfaces," IEEE Trans. on PAMI,
% 2012. (Accepted)

function newBSSFO = resamplingBSSFO( oldBSSFO )

accWeight = cumsum( oldBSSFO.weight );

for i=1:oldBSSFO.numBands
    oldBSSFO.selected( i ) = 0;
end

newBSSFO = oldBSSFO;
for i=1:oldBSSFO.numBands
    idx = find( accWeight > rand );
    newBSSFO.sample(:, i) = oldBSSFO.sample( :, idx(1) );
    if oldBSSFO.selected( idx(1) ) == 0
        oldBSSFO.selected( idx(1) ) = 1;
        newBSSFO.selected( i ) = 1;   % keep current state, do not add random noise
    end
end

newBSSFO.weight = ones(newBSSFO.numBands, 1)./ newBSSFO.numBands;

author: PatternRecognition

repo: OpenBMI


text:
function sF = times(sF1,sF2)

if isnumeric(sF1)
  sF = sF2;

  sF.values = sF2.values;

  sF.values = times(sF1, sF.values);

elseif isnumeric(sF2)

  sF = sF1;

  sF.values = sF1.values;

  sF.values = times(sF1.values, sF2);

else

  sF=sF1;

  sF.values = times(sF1.values, sF2.values);

end

end

author: mtex-toolbox

repo: mtex

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