text:
function fem2d_pack_test01 ( )

%*****************************************************************************80
%
%% TEST01 tests BANDWIDTH_MESH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  BANDWIDTH_MESH computes the geometric bandwidth:\n' );
  fprintf ( 1, '  of a finite element mesh.\n' );

  nelemx = 2;
  nelemy = 6;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  NELEMX = %d\n', nelemx );
  fprintf ( 1, '  NELEMY = %d\n', nelemy );

  element_order = 6;
  element_num = grid_element_num ( 'T6', nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  ELEMENT_ORDER = %d\n', element_order );
  fprintf ( 1, '  ELEMENT_NUM   = %d\n', element_num );

  element_node = grid_t6_element ( nelemx, nelemy );

  grid_print ( element_order, element_num, element_node );

  [ ml, mu, m ] = bandwidth_mesh ( element_order, element_num, ...
    element_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Lower bandwidth ML = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );
  fprintf ( 1, '  Total bandwidth M  = %d\n', m );

  return
end


text:
% ------------------------------------------------------------------------
%          Crop, Rotate, Adjust Contrast, and Downsample Histology
% ------------------------------------------------------------------------

%%  SET FILE AND PARAMETERS

% * remember to run one cell at a time, instead of the whole script at once *

% directory of histology images
image_folder = 'C:\Drive\Histology\cfos';

% directory to save the processed images -- can be the same as the above image_folder
% results will be put inside a new folder called 'processed' inside of this image_folder
save_folder = 'C:\Drive\Histology\cfos';

% name of images, in order anterior to posterior or vice versa
% once these are downsampled they will be named ['original name' '_processed.tif']
image_file_names = dir([image_folder filesep '*.tif']); % get the contents of the image_folder
image_file_names = natsortfiles({image_file_names.name});
% image_file_names = {'slide no 2_RGB.tif','slide no 3_RGB.tif','slide no 4_RGB.tif'}; % alternatively, list each image in order

% if the images are individual slices (as opposed to images of multiple
% slices, which must be cropped using the cell CROP AND SAVE SLICES)
image_files_are_individual_slices = true;

% use images that are already at reference atlas resolution (here, 10um/pixel)
use_already_downsampled_image = true;

% pixel size parameters: microns_per_pixel of large images in the image
% folder (if use_already_downsampled_images is set to false);
% microns_per_pixel_after_downsampling should typically be set to 10 to match the atlas
microns_per_pixel = 3.233;
microns_per_pixel_after_downsampling = 10;

% ----------------------
% additional parameters
% ----------------------

% if the images are cropped (image_file_are_individual_slices = false),
% name to save cropped slices as; e.g. the third cropped slice from the 2nd
% image containing many slices will be saved as: save_folder/processed/save_file_name02_003.tif
save_file_name = 'SS096_';

% increase gain if for some reason the images are not bright enough
gain = 1;

% plane to view ('coronal', 'sagittal', 'transverse')
plane = 'coronal';

% size in pixels of reference atlas brain. For coronal slice, this is 800 x 1140
if strcmp(plane,'coronal')
    atlas_reference_size = [800 1140];
elseif strcmp(plane,'sagittal')
    atlas_reference_size = [800 1320];
elseif strcmp(plane,'transverse')
    atlas_reference_size = [1140 1320];
end

% finds or creates a folder location for processed images --
% a folder within save_folder called processed
folder_processed_images = fullfile(save_folder, 'processed');
if ~exist(folder_processed_images)
    mkdir(folder_processed_images)
end

%% LOAD AND PROCESS SLICE PLATE IMAGES

% close all figures
close all

% if the images need to be downsampled to 10um pixels (use_already_downsampled_image = false),
% this will downsample and allow you to adjust contrast of each channel of each image from image_file_names
%
% if the images are already downsampled (use_already_downsampled_image = true), this will allow
% you to adjust the contrast of each channel
%
% Open Histology Viewer figure
try; figure(histology_figure);
catch; histology_figure = figure('Name','Histology Viewer'); end
warning('off', 'images:initSize:adjustingMag'); warning('off', 'MATLAB:colon:nonIntegerIndex');

% Function to downsample and adjust histology image
HistologyBrowser(histology_figure, save_folder, image_folder, image_file_names, folder_processed_images, image_files_are_individual_slices, ...
            use_already_downsampled_image, microns_per_pixel, microns_per_pixel_after_downsampling, gain)

%% CROP AND SAVE SLICES -- run once the above is done, if image_file_are_individual_slices = false

% close all figures
close all

% run this function if the images from image_file_names have several
% slices, per image (e.g. an image of an entire histology slide)
%
% this function allows you to crop all the slices you would like to
% process im each image, by drawing rectangles around them in the figure.
% these will be further processed in the next cell
if ~image_files_are_individual_slices
    histology_figure = figure('Name','Histology Viewer');
    HistologyCropper(histology_figure, save_folder, image_file_names, atlas_reference_size, save_file_name, use_already_downsampled_image)
else
    disp('individually cropped slices already available')
end

%% GO THROUGH TO FLIP HORIZONTAL SLICE ORIENTATION, ROTATE, SHARPEN, and CHANGE ORDER

% close all figures
close all

% this takes images from folder_processed_images ([save_folder/processed]),
% and allows you to rotate, flip, sharpen, crop, and switch their order, so they
% are in anterior->posterior or posterior->anterior order, and aesthetically pleasing
%
% it also pads images smaller than the reference_size and requests that you
% crop images larger than this size
%
% note -- presssing left or right arrow saves the modified image, so be
% sure to do this even after modifying the last slice in the folder
slice_figure = figure('Name','Slice Viewer');
SliceFlipper(slice_figure, folder_processed_images, atlas_reference_size)

author: cortex-lab

repo: allenCCF

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


text:
function showgraph(node,edge)
%% SHOWGRAPH displays a planar graph
%
%    showgraph(node,edge) displays a planar undirected graph.
%
%   See also showmesh, findedge
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

line([node(edge(:,1),1)'; node(edge(:,2),1)'],...
     [node(edge(:,1),2)'; node(edge(:,2),2)'],...
     'LineWidth',1,'Color',[0.125 0.5 0.125]);
hold on
plot(node(:,1),node(:,2),'k.', 'MarkerSize', 12);
axis equal; axis off

repo: Aero-matlab


text:
% Example usage of PoissonDistributionX

lambda = 5;

pd = PoissonDistributionX(lambda);

mu = pd.Mean;
covar = pd.Covar;
samples = pd.random(5000000);

lik = pd.pdf(samples);

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function lambda = carry_eigenvalues ( n, alpha )

%*****************************************************************************80
%
%% CARRY_EIGENVALUES returns the eigenvalues of the CARRY matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ALPHA, the numeric base being used in the addition.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  for i = 1 : n
    lambda(i,1) = 1.0 / alpha^(i-1);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%  quantization -- discretizes a set of random variates
%
%  ::
%
%    [x,w] = quantization(n,xw_discrete)
%
%  Args:
%
%     n (vector): vector containing the number of discrete quantities in
%        each dimension of the gaussian shocks
%
%     xw_discrete (1 x 2 cell array): The first cell is the vector of all
%        possible values of the discrete distribution. The second cell is
%        the vector of weights of each of those possible values
%
%  Returns:
%     :
%
%     - **x** [matrix]: grid combinations of the different variates, with
%        the number of shocks in rows and the number of combinations in
%        columns
%
%     - **w** [vector]: probability distribution over the different
%        combinations
%
%  Note:
%
%     - Only works with gaussian shocks at the moment.
%
%

author: jmaih


text:
function [a,n,r,epe,eph,epw] = mg_q1cd(xy,ev)
%mg_q1cd   convection-diffusion matrix generator for GMG
%   [a,n,r,epe,eph,epw] = mg_q1cd(xy,ev)
%   input
%           xy         vertex coordinate vector
%           ev         element mapping matrix
%   output
%           a          discrete diffusion operator
%           n          discrete convection operator
%           r          dummy variable
%           epe        viscosity normalised element peclet numbers
%           eph        flow specific element lengths
%           epw        centroid evaluated wind
%
%   IFISS function: DJS; 5 January 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

%   Analogous to femq1_cd
%
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));
%
% initialise global matrices
      a = sparse(nvtx,nvtx);
      n = sparse(nvtx,nvtx);
      r = sparse(nvtx,nvtx);
%
% 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);
	    ne = zeros(nel,4,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);
         [flowx,flowy] = gauss_transprt(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(:);
               ne(:,i,j) = ne(:,i,j) + flowx(:) .* phi(:,i) .* dphidx(:,j);
               ne(:,i,j) = ne(:,i,j) + flowy(:) .* phi(:,i) .* dphidy(:,j);
               end
	    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);
          n = n + sparse(nrow,ncol,ne(:,krow,kcol),nvtx,nvtx);
          end
      end
%
%
% computation of element Peclet number (at the centroid)
% rectangle specific calculation here
      hx=abs(xl_v(:,2)-xl_v(:,1)); hy=abs(yl_v(:,3)-yl_v(:,2));
      [flowx,flowy] = gauss_transprt(0,0,xl_v,yl_v);
      flow_l2 = sqrt(flowx(:) .* flowx(:) + flowy(:) .* flowy(:));
        if all(flowx==0), flow_h=hy;
	elseif all(flowy==0), flow_h=hx;
		else
          angle = atan(abs(flowy./flowx));
          flow_h = min([hx./cos(angle),hy./sin(angle)],[],2);
        end
	  eph = flow_h;
      epe = flow_h.*flow_l2/2;
	  epw = flow_l2;
%
return


text:
function [x, y] = CheckRC(x, y, im)
y = max(y, 1);
y = min(y, size(im, 1));
x = max(x, 1);
x = min(x, size(im, 2));

author: decouples


text:
function btv_test05 ( )

%*****************************************************************************80
%
%% BTV_TEST05 tests BURGERS_TIME_VISCOUS with the expansion initial condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTV_TEST05\n' );
  fprintf ( 1, '  Test BURGERS_TIME_VISCOUS with the expansion initial condition.\n' );
  fprintf ( 1, '  Use periodic boundaries.\n' );

  nx = 81;
  nt = 200;
  t_max = 2.0;
  nu = 0.01;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial condition: expansion\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Viscosity = %g\n', nu );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_viscous ( @ic_expansion, nx, nt, t_max, nu, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 5 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers equation solutions over time, initial condition expansion' )

  filename = 'btv_test05.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end


text:
function title = p08_title ( )

%*****************************************************************************80
%
%% P08_TITLE returns the title of problem 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Extrapolation test data';

  return
end


text:
function [adfreq, n, ts, fn] = plx_ad_gap_info(filename, channel)
% plx_ad_gap_info(filename, channel): read a/d info from a .plx or .pl2 file
%
% [adfreq, n, ts, fn] = plx_ad_gap_info(filename, channel)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   channel - 0-based channel number or channel name
%
%           a/d data come in fragments. Each fragment has a timestamp
%           and a number of a/d data points. The timestamp corresponds to
%           the time of recording of the first a/d value in this fragment.
%           All the data values stored in the vector ad.
%
% OUTPUT:
%   adfreq - digitization frequency for this channel
%   n - total number of data points
%   ts - array of fragment timestamps (one timestamp per fragment, in seconds)
%   fn - number of data points in each fragment

adfreq = 0;
n = 0;
ts = -1;
fn = -1;

if nargin ~= 2
    error 'expected 2 input arguments';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    pl2Channel = channel;
    % PL2 uses 1-based channel numbers
    if ischar(channel) == 0
        pl2Channel = channel+1;
    end
    pl2ad = PL2Ad(filename, pl2Channel);
    if numel(pl2ad.Values) > 0
        adfreq = pl2ad.ADFreq;
        n = numel(pl2ad.Values);
        ts = pl2ad.FragTs;
        fn = pl2ad.FragCounts;
    end
    return;
end

channelNumber = plx_ad_resolve_channel(filename, channel);
if channelNumber == -1
    fprintf('\n plx_ad_gap_info: no header for the specified A/D channel.');
    fprintf('\n    use plx_ad_info(filename) to print the list of a/d channels\n');
    return
end

[adfreq, n, ts, fn] = mexPlex(24, filename, channelNumber);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function test_suite = test_fmri_deoblique
% tests for cosmo_fmri_deoblique
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_fmri_deoblique_basics
    ds=cosmo_synthetic_dataset('size','normal','ntargets',1,'nchunks',1);
    % make dataset oblique (manually)
    ds.a.vol.mat(1,1)=.8;
    ds.a.vol.mat(2,1)=.6;

    ds_deoblique=cosmo_fmri_deoblique(ds);

    mat=eye(4);
    mat(2,2)=2;
    mat(3,3)=2;
    mat(1:3,4)=[-3.2 -2.4 -3];

    assertEqual(ds_deoblique.a.vol.mat,mat);
    ds.a.vol.mat=mat;
    assertEqual(ds,ds_deoblique);

    ds_deoblique2=cosmo_fmri_deoblique(ds_deoblique);
    assertEqual(ds_deoblique,ds_deoblique2);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% DO VARIOUS ROBOTS: 2 LINK, HYPER2D, HYPER3D

% Inverse kinematics is the problem of finding the robot joint coordinates,
% given a homogeneous transform representing the last link of the manipulator.
% It is very useful when the path is planned in Cartesian space, for instance
% a straight line path as shown in the trajectory demonstration.
%
% First generate the transform corresponding to a particular joint coordinate,

mdl_puma560
q = [0 -pi/4 -pi/4 0 pi/8 0]
T = p560.fkine(q)

% Now the inverse kinematic procedure for any specific robot can be derived
% symbolically and in general an efficient closed-form solution can be
% obtained.  However we are given only a generalized description of the
% manipulator in terms of kinematic parameters so an iterative solution will
% be used. The procedure is slow, and the choice of starting value affects
% search time and the solution found, since in general a manipulator may
% have several poses which result in the same transform for the last
% link. The starting point for the first point may be specified, or else it
% defaults to zero (which is not a particularly good choice in this case)

qi = p560.ikine(T);

% and in fact it does not converge

qi

% We can help the solution along by using the 'pinv' option

qi = p560.ikine(T, 'pinv');
% and the result
qi
% is the same as the original set of joint angles
q

% However in general this will not be the case, there are multiple
% solutions, and the solution that is found depends on the initial
% choice of angles.
%
% A more efficient approach is to use an analytic solution and the toolbox
% supports the common case of a 6-axis robot arm with a spherical wrist

qi = p560.ikine6s(T)
% which is different to the original joint angles, but as expected

p560.fkine(qi)

% it does give the same end-effector pose.

% The analytic solution allows the specific solution to be specified
% using a character string and to get the same set of joint angles

p560.ikine6s(T, 'rdf')
% where we have specified that the robot is in a right-handed configuration
% (r), with its elbow down (d), and the wrist flipped (f).

% A solution is not always possible, for instance if the specified
% transform describes a point out of reach of the manipulator.  As
% mentioned above the solutions are not necessarily unique, and there
% are singularities at which the manipulator loses degrees of freedom
% and joint coordinates become linearly dependent.

% Inverse kinematics may also be computed for a trajectory.
% If we take a Cartesian straight line path between two poses in 50 steps

T1 = transl(0.6, -0.5, 0.0) % define the start point
T2 = transl(0.4, 0.5, 0.2)	% and destination
T = ctraj(T1, T2, 50); 	% compute a Cartesian path

% now solve the inverse kinematics

q = p560.ikine6s(T);
about q
% which has one row per time step and one column per joint angle

% Let's examine the joint space trajectory that results in straightline
% Cartesian motion

subplot(3,1,1); plot(q(:,1)); xlabel('Time (s)'); ylabel('Joint 1 (rad)');
subplot(3,1,2); plot(q(:,2)); xlabel('Time (s)'); ylabel('Joint 2 (rad)');
subplot(3,1,3); plot(q(:,3)); xlabel('Time (s)'); ylabel('Joint 3 (rad)');

% This joint space trajectory can now be animated
clf
p560.plot(q)

author: petercorke

repo: robotics-toolbox-matlab


y = 17;
x = (y + 3) / 2);


text:
function z=mergeSortedArrays(x,y)
%MERGESORTEDARRAYS Given two arrays that are each sorted in ascending
%                  order, merge them into one large array that is in
%                  ascending order.
%
%INPUTS: x An xDimX1 or 1XxDim array that has been sorted in ascending
%          order.
%        y A yDimX1 or 1XyDim array that has been sorted in ascending
%          order.
%
%OUTPUTS: z An (xDim+yDim)X1 array containing all of the elements of x and
%           y, in ascending order.
%
%The algorthm is taken from Chapter 5.2.4 of [1].
%
%REFERENCES:
%[1] D. Knuth, The Art of Computer Programming: Sorting and Searching, 2nd
%    ed. Reading, MA: Addison-Wesley, 1998, vol. 3.
%
%December 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

m=length(x);
n=length(y);
zLen=n+m;
z=zeros(zLen,1);

%Initialize
i=1;
j=1;
k=1;

%Find smaller
while(1)
    if(x(i)<=y(j))
        z(k)=x(i);
        k=k+1;
        i=i+1;
        if(i<=m)
            continue;
        else
            z(k:zLen)=y(j:n);
            break;
        end
    else
        z(k)=y(j);
        k=k+1;
        j=j+1;

        if(j<=n)
            continue;
        else
           z(k:zLen)=x(i:m);
           break;
        end
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function f = div( f )
%DIV   Numerical divergence of a DISKFUNV.
%   D = DIVERGENCE( F ) returns the numerical divergence of the
%   DISKFUNV.
%
%  This is shorthand for the command DIVERGENCE.
%
% See also DIVERGENCE, GRAD, CURL.

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

f = divergence( f );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function newlp = mrRefreshInplanes(lp,obXM,obYM,curInplane,deleteFlag)
% function newlp = mrRefreshInplanes(lp,obXM,obYM,curInplane,deleteFlag)
% ---------------------------------------------------
% Function: refreshes the inplane grid image replacing old handles with
% new ones derived from obXM and obYM
% Author: SPG
% Date: 2.26.97
% Notes: object handling in matlab is seriously screwed up. So be wary of
% altering any delete(lp) calls here or elsewhere, because I don't think
% logic applies to them.

global sagwin

figure(sagwin);

if deleteFlag == 1
  for i=1:length(lp)
    delete(lp(i));
  end
end

if ~isempty(lp)
 for i=1:length(lp)
   if i == curInplane
        newlp(i)=line(obXM(i,:),obYM(i,:),'Color','r');
   else
        newlp(i)=line(obXM(i,:),obYM(i,:),'Color','b');
   end
  end
else
 newlp = [];
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function f2 = p06_f2 ( x )

%*****************************************************************************80
%
%% P06_F2 evaluates the second derivative for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    LE Scales,
%    Introduction to Non-Linear Optimization,
%    Springer, 1985.
%
%  Parameters:
%
%    Input, real X, the values of the variables.
%
%    Output, real F2, the second derivative.
%
  f2 = 0.0;

  return
end


text:
function biosigpathfirst

str2doublepath = fileparts( which('str2double') );
sopenpath      = fileparts( which('sopen') );
if ~strcmp(str2doublepath,sopenpath)
    rmpath(sopenpath);
    addpath(sopenpath,'-begin');
end


text:
function res = ctranspose(obj,useconj);
    if nargin == 1
        useconj = 1;
    end

    if length(obj.dimensions) > 2
        error('Cannot transpose array');
    end

    % make new memory mapped data file
    % --------------------------------
    newFileName = mmo.getnewfilename;
    copyfile(obj.dataFile, newFileName);

    res = obj;
    res.dimensions = [ obj.dimensions(2) obj.dimensions(1) ];
    res.dataFile = newFileName;
    tmpMMO1 = memmapfile(obj.dataFile, 'writable', obj.writable, 'format', { 'single' obj.dimensions 'x' });
    tmpMMO2 = memmapfile(res.dataFile, 'writable',         true, 'format', { 'single' res.dimensions 'x' });

    % copy the data
    % -------------
    if length(obj.dimensions) == 1 || obj.dimensions(1) > obj.dimensions(2)
        for index = 1:size(obj,2)
            s.type = '()';
            s.subs = { ':' index };
            if useconj, tmpMMO2.Data.x(index,:) = conj(subsref(tmpMMO1.Data.x,s));
            else        tmpMMO2.Data.x(:,index) =      subsref(tmpMMO1.Data.x,s);
            end
        end
    else
        for index = 1:size(obj,1)
            s.type = '()';
            s.subs = { index ':' };
            if useconj, tmpMMO2.Data.x(:,index) = conj(subsref(tmpMMO1.Data.x,s));
            else        tmpMMO2.Data.x(:,index) =      subsref(tmpMMO1.Data.x,s);
            end
        end
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [x,b,a] = bst_bandpass_filtfilt(x, Fs, HighPass, LowPass, isStopBand, FilterType)
% BST_BANDPASS_FILTFILT: Bandpass filter for the signal x, using the filtfilt function (used by default after Nov 2014)
%
% USAGE:  [x,b,a] = bst_bandpass_filtfilt(x, Fs, HighPass, LowPass, isStopBand=0, FilterType='fir')
%
% INPUT:
%    - x          : [nChannels,nTime] signal to process
%    - Fs         : Sampling frequency
%    - HighPass   : Frequency below this value are filtered (set to 0 for low-pass filter only)
%    - LowPass    : Frequency above this value are filtered (set to 0 for high-pass filter only)
%    - isStopBand : If 1, create a stop-band filter instead of a pass-band filter
%    - FilterType : 'fir' or 'iir'
%
% OUTPUT:
%    - x   : Filtered signals
%    - b,a : Filter coefficients, as defined in all the Matlab functions
%
% DESCRIPTION:
%    - A linear phase FIR filter is created and applied both forward and backward using filtfilt.
%    - Function "filtfilt" is used to employ the filtering "mirror" trick, and to reset automatically the group delay.
%    - Function "kaiserord" and "kaiser" are used to set the necessary order for fir1.
%    - The transition band is hard-coded.
%    - Requires Signal Processing Toolbox for the following functions: kaiserord, kaiser, ellipord,
%      If not, using Octave-based alternatives

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

% ===== PARSE INPUTS =====
if (nargin < 6) || isempty(FilterType)
    FilterType = 'fir';
end
if (nargin < 5) || isempty(isStopBand)
    isStopBand = 0;
end
if isempty(HighPass)
    HighPass = 0;
end
if isempty(LowPass)
    LowPass = 0;
end
% If both high-pass and low-pass are zero: return signals unaltered
if (HighPass == 0) && (LowPass == 0)
    disp('BST_BANDPASS> Error: No frequency band in input');
    return;
end

% ===== FILTER PARAMETERS =====
PASSBAND_RIPPLE = 5;    % percent pass band ripple
PASSBAND_DB     = 1;    % dB of pass band ripple
STOP_ATTEN_DB   = 40;   % dB of attenuation in the stop band
TRANSITION_BAND = 0.05; % normalized to Nyquist, the allowed transition band
% We use filtfilt, which doubles the effect (attenuation and ripple)
PASSBAND_RIPPLE = PASSBAND_RIPPLE/2;
STOP_ATTEN_DB   = STOP_ATTEN_DB/2;
% Conversion from percent
Ripple = PASSBAND_RIPPLE/100;
% Stop band attenuation
Atten  = 10^(-STOP_ATTEN_DB/20);

% Convert frequencies to normalized form
Nyquist = Fs/2;
f_highpass = HighPass / Nyquist;
f_lowpass  = LowPass  / Nyquist;
% Reasonable digital transition band
f_highstop = f_highpass - TRANSITION_BAND;
f_lowstop  = f_lowpass  + TRANSITION_BAND;

% ===== CREATE FILTER =====
switch FilterType
    % ===== FIR =====
    case 'fir'
        % Build the general case first
        fcuts = [f_highstop, f_highpass, f_lowpass, f_lowstop];
        % Stop-band
        if isStopBand
            mags = [1 0 1];               % filter magnitudes
            devs = [Ripple Atten Ripple]; % deviations
        % Pass-band
        else
            mags = [0 1 0];               % filter magnitudes
            devs = [Atten Ripple Atten];  % deviations
        end
        % Now adjust for desired properties
        fcuts = max(0,fcuts);     % Can't go below zero
        fcuts = min(1-eps,fcuts); % Can't go above or equal to 1
        % We have implicitly created a bandpass, but now adjust for desired filter
        if (f_lowpass == 0)  % User didn't want a lowpass
            fcuts(3:4) = [];
            mags(3) = [];
            devs(3) = [];
        end
        if (f_highpass == 0)  % User didn't want a highpass
            fcuts(1:2) = [];
            mags(1) = [];
            devs(1) = [];
        end
        % Generate FIR filter
        [n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,2);
        n = n + rem(n,2);  % ensure even order
        b = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
        a = 1;

    % ===== IIR =====
    case 'iir'
        % Stop-band
        if isStopBand
            ftype = 'stop';
            Ws = [f_highpass f_lowpass]; % the range of stopped
            Wp = [f_highstop f_lowstop]; % the transition band
        % Pass-band
        else
            ftype = 'bandpass';
            Ws = [f_highstop f_lowstop]; % the transition band
            Wp = [f_highpass f_lowpass]; % the passband
        end
        % Now handle extremes
        Ws = max(eps,Ws);   % Can't be zero or less
        Wp = max(eps,Wp);
        Ws = min(1-eps,Ws); % Can't be 1
        Wp = min(1-eps,Wp);
        % Now handle highpass or lowpass only
        if (f_lowpass == 0)  % User didn't want a lowpass
            ftype = 'high';
            Ws(2) = [];
            Wp(2) = [];
        end
        if (f_highpass == 0)  % User didn't want a highpass
            ftype = 'low';
            Ws(1) = [];
            Wp(1) = [];
        end
        % Generate IIR filter
        [n,WP] = ellipord(Wp,Ws,PASSBAND_DB,STOP_ATTEN_DB);
        [b,a]  = ellip(n,PASSBAND_DB,STOP_ATTEN_DB,WP,ftype);
end

% ===== FILTER THE DATA =====
Ntime = size(x,2);
% Remove the mean of the data before filtering, which wrecks most filters
xmean = mean(x,2);
x = bst_bsxfun(@minus, x, xmean)';    % Transposed output (time is now down the columns)
% Using filtfilt to use the mirroring trick and remove group delay
try
    x = filtfilt(b,a,x)';     % Transposed output

    % OCTAVE IMPLEMENTATION
    % http://octave-signal.sourcearchive.com/documentation/1.0.8/filtfilt_8m-source.html
catch
    fprintf('Sequence too short for filtfilt, using alternate approach.\n')
    xmirror = [flipud(x); x; flipud(x)];      % Mirror either end
    xmirror = filter(b,1,xmirror);            % Filter
    xmirror = flipud(xmirror);                % Reverse in time
    xmirror = flipud(filter(b,1,xmirror));    % Filter and flip again
    x = xmirror(Ntime + (1:Ntime),:)';        % Transposed output
end
% Restore the mean of the signal (only if there is no high-pass filter)
if (f_highpass == 0)
    x = bst_bsxfun(@plus, x, xmean);
end

author: brainstorm-tools


text:
function [X,Y,vals,labI]=mp_conic(optn,varargin)
% MP_CONIC  Conic projections
%           This function should not be used directly; instead it is
%           is accessed by various high-level functions named M_*.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 2/Apr/1997
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
% Mathematical formulas for the projections and their inverses are taken from
%
%      Snyder, John P., Map Projections used by the US Geological Survey,
%      Geol. Surv. Bull. 1532, 2nd Edition, USGPO, Washington D.C., 1983.
%
% These are  conic projections with two standard parallels, useful
% for showing limited areas at mid-latitudes.
%    Albers equal-area - has an equal-area property
%    Lambert conformal - is conformal
%
%  7/6/99 - fixed tendency to re-define .ulongs if .clong set by user
%  3/4/02 - added error if parallels are equidistant from equator (i.e. not conic projection really)
%  24/10/08 - added ellipsoidal earth computations for lambert conformal conic projection
%  16/10/09 - added ellipsoidal earth computations for albers conic projection
%  06/03/17 - changed 'false_origin' to 'origin' as option name, also it
%             didn't work correctly with the 'normal' spheroid so this was
%             fixed, AND changed the default parallels so they were at 25%
%             and 75% limits instead of being a single parallel at the
%             center to prevent blowups with the albers ellipsoidal
%             projection AND made it work with SPHERE ellipsoid.

global MAP_PROJECTION MAP_VAR_LIST

MAP_ELLIP=mc_ellips;

name={'Albers Equal-Area Conic','Lambert Conformal Conic'};

pi180=pi/180;

switch optn

  case 'name'

     X=name;

  case {'usage','set'}

     m_names=fieldnames(MAP_ELLIP);

     X=char({['     ''' varargin{1} ''''],...
              '     <,''lon<gitude>'',[min max]>',...
              '     <,''lat<itude>'',[min max]>',...
              '     <,''clo<ngitude>'',value>',...
              '     <,''par<allels>'',[lat1 lat2]>',...
              '     <,''rec<tbox>'', ( ''on'' | ''off'' )>',...
              '     <,''ell<ipsoid>'', one of',...
    reshape(sprintf('         %6s',m_names{:}),15,length(m_names))',...
        '                               >',...
	      '     <,''ori<gin>'', [long lat]>'});

  case 'get'

     X=char([' Projection: ' MAP_PROJECTION.name '  (function: ' MAP_PROJECTION.routine ')'],...
            [' longitudes: ' num2str(MAP_VAR_LIST.ulongs) ' (centered at ' num2str(MAP_VAR_LIST.clong) ')'],...
            [' latitudes: ' num2str(MAP_VAR_LIST.ulats) ],...
            [' standard parallels: ' num2str(MAP_VAR_LIST.parallels) ],...
            [' Rectangular border: ' MAP_VAR_LIST.rectbox ],...
            [' ellipsoid: ' MAP_VAR_LIST.ellipsoid ],...
            [' origin: ' num2str(MAP_VAR_LIST.origin) ]);

  case 'initialize'

    MAP_VAR_LIST=[];
    MAP_PROJECTION.name=varargin{1};
    MAP_VAR_LIST.ulongs=[-180 -50];
    MAP_VAR_LIST.ulats=[10 85];
    MAP_VAR_LIST.parallels=NaN;
    MAP_VAR_LIST.clong=NaN;
    MAP_VAR_LIST.origin=NaN;
    MAP_VAR_LIST.rectbox='off';
    MAP_VAR_LIST.ellipsoid = 'normal';
    MAP_VAR_LIST.aussiemode=false;
    k=2;longs_def=0;
    while k<length(varargin)
      switch varargin{k}(1:3)
         case 'lon'
           MAP_VAR_LIST.ulongs=varargin{k+1}(:)';longs_def=1;
           if MAP_VAR_LIST.ulongs(1)>MAP_VAR_LIST.ulongs(2)
              MAP_VAR_LIST.ulongs=MAP_VAR_LIST.ulongs([2 1]);
           end
         case 'clo'
           MAP_VAR_LIST.clong=varargin{k+1};
         case 'lat'
           MAP_VAR_LIST.ulats=varargin{k+1}(:)';
         case 'par'
           MAP_VAR_LIST.parallels=varargin{k+1};
         case 'rec'
             switch lower(varargin{k+1}(1:2))
                 case {'on','bo'}
                    MAP_VAR_LIST.rectbox='on';
                 case 'of'
                    MAP_VAR_LIST.rectbox='on';
                 otherwise
                     error(['m_proj: Unrecognized box option: ' varargin{k+1}]);
             end
         case 'ell'
           MAP_VAR_LIST.ellipsoid=varargin{k+1};
         case 'ori'
           MAP_VAR_LIST.origin=varargin{k+1};
         case 'fal'
           error(' FALSE_ORIGIN option has been renamed to ORIGIN - Change your m_proj call! ');
	 case 'aus'  % aussiemode - my joke
	   if strcmp(varargin{k+1},'on')
	     MAP_VAR_LIST.aussiemode=true;
	   end
         otherwise
           disp(['Unknown option: ' varargin{k}]);
      end
       k=k+2;
    end
    if isnan(MAP_VAR_LIST.clong)
      if isnan(MAP_VAR_LIST.origin)
         MAP_VAR_LIST.clong=mean(MAP_VAR_LIST.ulongs);
      else
         MAP_VAR_LIST.clong=MAP_VAR_LIST.origin(1);
      end
    elseif  ~longs_def
        MAP_VAR_LIST.ulongs=MAP_VAR_LIST.clong+[-180 180];
    end
    if isnan(MAP_VAR_LIST.parallels), MAP_VAR_LIST.parallels=mean(MAP_VAR_LIST.ulats)*[1 1]+diff(MAP_VAR_LIST.ulats)*[-1/6 1/6]; end % change default mar/2017
    if isnan(MAP_VAR_LIST.origin), MAP_VAR_LIST.origin=[MAP_VAR_LIST.clong mean(MAP_VAR_LIST.parallels)]; end

    MAP_VAR_LIST.rlongs=MAP_VAR_LIST.ulongs*pi180;
    MAP_VAR_LIST.rlats=MAP_VAR_LIST.ulats*pi180;
    MAP_VAR_LIST.rparallels=MAP_VAR_LIST.parallels*pi180;
    MAP_VAR_LIST.rorigin=MAP_VAR_LIST.origin*pi180;

    MAP_VAR_LIST.ellip=getfield(MAP_ELLIP,MAP_VAR_LIST.ellipsoid);

    % These are constants used by the projection formulas

    switch MAP_PROJECTION.name
      case name(1)
        if MAP_VAR_LIST.ellip(2)==0 % spherical
           MAP_VAR_LIST.n=sum(sin(MAP_VAR_LIST.rparallels))/2;
           if MAP_VAR_LIST.n==0
               error('Your parallels are equidistant from the equator - use a cylindrical projection!');
           end
           MAP_VAR_LIST.C=cos(MAP_VAR_LIST.rparallels(1)).^2+2*MAP_VAR_LIST.n*sin(MAP_VAR_LIST.rparallels(1));
           MAP_VAR_LIST.rho0=MAP_VAR_LIST.ellip(1)*sqrt(MAP_VAR_LIST.C-2*MAP_VAR_LIST.n*sin( MAP_VAR_LIST.rorigin(2)   ))/MAP_VAR_LIST.n;
        else
 	       e=sqrt(2*MAP_VAR_LIST.ellip(2)-MAP_VAR_LIST.ellip(2)^2);
	       m12=cos(MAP_VAR_LIST.rparallels)./sqrt(1-(e.*sin(MAP_VAR_LIST.rparallels)).^2);
	       q12=(1-e.^2)*(sin(MAP_VAR_LIST.rparallels)./(1-(e.*sin(MAP_VAR_LIST.rparallels)).^2) - ...
	            1./(2*e)*log((1-e.*sin(MAP_VAR_LIST.rparallels))./(1+e.*sin(MAP_VAR_LIST.rparallels))) );
	       q0= (1-e.^2)*(sin(MAP_VAR_LIST.rorigin(2))./(1-(e.*sin(MAP_VAR_LIST.rorigin(2))).^2) - ...
	            1./(2*e)*log((1-e.*sin(MAP_VAR_LIST.rorigin(2)))./(1+e.*sin(MAP_VAR_LIST.rorigin(2)))) );
           if diff( MAP_VAR_LIST.rparallels )==0
	          MAP_VAR_LIST.n=sin(MAP_VAR_LIST.rparallels(1));
           else
	          MAP_VAR_LIST.n=-diff(m12.^2)/diff(q12);
           end
	       MAP_VAR_LIST.C=m12(1).^2 + MAP_VAR_LIST.n.*q12(1);
	       MAP_VAR_LIST.rho0=MAP_VAR_LIST.ellip(1)*sqrt(MAP_VAR_LIST.C-MAP_VAR_LIST.n*q0)/MAP_VAR_LIST.n;
        end
      case name(2)
        if strcmp(MAP_VAR_LIST.ellipsoid,'normal')
           if diff(MAP_VAR_LIST.parallels)==0
             MAP_VAR_LIST.n=sin(MAP_VAR_LIST.rparallels(1));
           else
             MAP_VAR_LIST.n=-diff(log(cos(MAP_VAR_LIST.rparallels)))/diff(log(tan(MAP_VAR_LIST.rparallels/2+pi/4)));
           end
           MAP_VAR_LIST.F=cos(MAP_VAR_LIST.rparallels(1))/MAP_VAR_LIST.n* ...
                	  tan(pi/4+MAP_VAR_LIST.rparallels(1)/2).^MAP_VAR_LIST.n;
           MAP_VAR_LIST.rho0=MAP_VAR_LIST.F/tan(pi/4+ MAP_VAR_LIST.rorigin(2)/2).^MAP_VAR_LIST.n;
        else
 	      e=sqrt(2*MAP_VAR_LIST.ellip(2)-MAP_VAR_LIST.ellip(2)^2);
	      m12=cos(MAP_VAR_LIST.rparallels)./sqrt(1-(e.*sin(MAP_VAR_LIST.rparallels)).^2);
	      t12=tan(pi/4-MAP_VAR_LIST.rparallels/2)./( (1-e*sin(MAP_VAR_LIST.rparallels))./(1+e*sin(MAP_VAR_LIST.rparallels)) ).^(e/2);
	      tF=tan(pi/4-MAP_VAR_LIST.rorigin(2)/2)./( (1-e*sin(MAP_VAR_LIST.rorigin(2)))./(1+e*sin(MAP_VAR_LIST.rorigin(2))) ).^(e/2);
	      if diff(MAP_VAR_LIST.rparallels)==0
	         MAP_VAR_LIST.n=sin(MAP_VAR_LIST.rparallels(1));
          else
             MAP_VAR_LIST.n=diff(log(m12))/diff(log(t12));
          end
	      MAP_VAR_LIST.F=m12(1)/MAP_VAR_LIST.n/t12(1).^MAP_VAR_LIST.n;
	      MAP_VAR_LIST.rho0=MAP_VAR_LIST.ellip(1)*MAP_VAR_LIST.F*tF.^MAP_VAR_LIST.n;
        end
    end

    % check for a valid ellipsoid. if not, use the normalized sphere

    if ~isfield(MAP_ELLIP,MAP_VAR_LIST.ellipsoid)
       MAP_VAR_LIST.ellipsoid = 'normal';
    end

    % Get X/Y and (if we are in a box) update the lat/long limits.

    mu_util('xylimits');
    if strcmp(MAP_VAR_LIST.rectbox,'on'),  mu_util('lllimits'); end

  case 'll2xy'

    long=varargin{1};
    lat=varargin{2};
    vals=zeros(size(long));

    % Clip out-of-range values (lat/long box)

    if ~strcmp(MAP_VAR_LIST.rectbox,'on') && ~strcmp(varargin{4},'off')
        vals=vals | long<=MAP_VAR_LIST.longs(1)+eps*10 | long>=MAP_VAR_LIST.longs(2)-eps*10 | ...
	              lat<=MAP_VAR_LIST.lats(1)+eps*10 |   lat>=MAP_VAR_LIST.lats(2)-eps*10;
        [long,lat]=mu_util('clip',varargin{4},long,MAP_VAR_LIST.longs(1),long<MAP_VAR_LIST.longs(1),lat);
        [long,lat]=mu_util('clip',varargin{4},long,MAP_VAR_LIST.longs(2),long>MAP_VAR_LIST.longs(2),lat);
        [lat,long]=mu_util('clip',varargin{4},lat,MAP_VAR_LIST.lats(1),lat<MAP_VAR_LIST.lats(1),long);
        [lat,long]=mu_util('clip',varargin{4},lat,MAP_VAR_LIST.lats(2),lat>MAP_VAR_LIST.lats(2),long);
    end

    switch MAP_PROJECTION.name
      case name(1)
        if MAP_VAR_LIST.ellip(2)==0   % spherical
           rho=MAP_VAR_LIST.ellip(1)*sqrt(MAP_VAR_LIST.C-2*MAP_VAR_LIST.n*sin(lat*pi180))/MAP_VAR_LIST.n;
        else
	       e=sqrt(2*MAP_VAR_LIST.ellip(2)-MAP_VAR_LIST.ellip(2)^2);
	       q= (1-e.^2)*(sin(lat*pi180)./(1-(e.*sin(lat*pi180)).^2) - ...
	            1./(2*e)*log((1-e.*sin(lat*pi180))./(1+e.*sin(lat*pi180))) );
	       rho=MAP_VAR_LIST.ellip(1)*sqrt(MAP_VAR_LIST.C-MAP_VAR_LIST.n*q)/MAP_VAR_LIST.n;
        end
      case name(2)
        if strcmp(MAP_VAR_LIST.ellipsoid,'normal')
           lat(lat==-90)=-89.999; % Prevents /0 problems in next line
           rho=MAP_VAR_LIST.F ./ tan(pi/4+lat*pi180/2).^MAP_VAR_LIST.n;
        else
	       e=sqrt(2*MAP_VAR_LIST.ellip(2)-MAP_VAR_LIST.ellip(2)^2);
	       t=tan(pi/4-lat*pi180/2)./( (1-e*sin(lat*pi180))./(1+e*sin(lat*pi180)) ).^(e/2);
	       rho=MAP_VAR_LIST.ellip(1)*MAP_VAR_LIST.F*t.^MAP_VAR_LIST.n;
        end
    end
    theta=MAP_VAR_LIST.n*(long-MAP_VAR_LIST.origin(1))*pi180;

    X=real(rho.*sin(theta));
    Y=real(MAP_VAR_LIST.rho0-rho.*cos(theta));

    % Clip out-of-range values (rectangular box)

    if strcmp(MAP_VAR_LIST.rectbox,'on') && ~strcmp(varargin{4},'off')
        vals= vals | X<=MAP_VAR_LIST.xlims(1)+eps*10 | X>=MAP_VAR_LIST.xlims(2)-eps*10 | ...
                     Y<=MAP_VAR_LIST.ylims(1)+eps*10 | Y>=MAP_VAR_LIST.ylims(2)-eps*10;
        [X,Y]=mu_util('clip',varargin{4},X,MAP_VAR_LIST.xlims(1),X<MAP_VAR_LIST.xlims(1),Y);
        [X,Y]=mu_util('clip',varargin{4},X,MAP_VAR_LIST.xlims(2),X>MAP_VAR_LIST.xlims(2),Y);
        [Y,X]=mu_util('clip',varargin{4},Y,MAP_VAR_LIST.ylims(1),Y<MAP_VAR_LIST.ylims(1),X);
        [Y,X]=mu_util('clip',varargin{4},Y,MAP_VAR_LIST.ylims(2),Y>MAP_VAR_LIST.ylims(2),X);
    end
    if MAP_VAR_LIST.aussiemode, Y=-Y; X=-X; end

  case 'xy2ll'

    pi180=pi/180;

    if MAP_VAR_LIST.aussiemode, varargin{2}=-varargin{2}; varargin{1}=-varargin{1}; end
    switch MAP_PROJECTION.name
      case name(1)
        rho=sqrt(varargin{1}.^2+(MAP_VAR_LIST.rho0-varargin{2}).^2);
        theta=atan(varargin{1}./(MAP_VAR_LIST.rho0-varargin{2}));
        if  MAP_VAR_LIST.ellip(2)==0
           Y=asin((MAP_VAR_LIST.C-(rho*MAP_VAR_LIST.n/MAP_VAR_LIST.ellip(1)).^2)/(2*MAP_VAR_LIST.n))/pi180;
        else
	       e=sqrt(2*MAP_VAR_LIST.ellip(2)-MAP_VAR_LIST.ellip(2)^2);
           q=(MAP_VAR_LIST.C - (rho.*MAP_VAR_LIST.n/MAP_VAR_LIST.ellip(1)).^2)./MAP_VAR_LIST.n;
	       % Y is computed iteratively
	       Y=asin(q/2);
	       for k=1:4
	          Y=Y+(1-(e.*sin(Y)).^2).^2./(2*cos(Y)).*( q./(1-e.^2) - sin(Y)./(1-(e.*sin(Y)).^2) + ...
	             1./(2*e)*log( (1-e*sin(Y))./(1+e*sin(Y)) ) );
           end
	       Y=Y/pi180;
        end
        % The pole is an arc in this projection, so for points inside that
        % arc there is no inverse. If so, the math above can return complex
        % values - instead set these to NaN
        if ~isreal(Y)
           Y(imag(Y)>1e-7)=NaN;
           Y=real(Y);
        end

      case name(2)
        rho=sign(MAP_VAR_LIST.n)*sqrt(varargin{1}.^2+(MAP_VAR_LIST.rho0-varargin{2}).^2);
        theta=atan(varargin{1}./(MAP_VAR_LIST.rho0-varargin{2}));
	    if strcmp(MAP_VAR_LIST.ellipsoid,'normal')
           Y=(2*atan((MAP_VAR_LIST.F./rho).^(1/MAP_VAR_LIST.n))-pi/2)/pi180;
        else
	       e=sqrt(2*MAP_VAR_LIST.ellip(2)-MAP_VAR_LIST.ellip(2)^2);
	       tp=(rho./MAP_VAR_LIST.ellip(1)/MAP_VAR_LIST.F).^(1./MAP_VAR_LIST.n);
	       % Y is computed iteratively
	       Y=pi/2 - 2*atan(tp);
	       for k=1:4
	         Y=pi/2 - 2*atan(tp.*((1-e*sin(Y))./(1+e*sin(Y))).^(e/2) );
           end
	       Y=Y/pi180;
        end
    end
    % Clip out-of-range values (lat/long box)
    X=MAP_VAR_LIST.origin(1)+(theta/MAP_VAR_LIST.n)/pi180;

  case 'xgrid'

    [X,Y,vals,labI]=mu_util('xgrid',MAP_VAR_LIST.longs,MAP_VAR_LIST.lats,varargin{1},3,varargin{2:3});

  case 'ygrid'

    [X,Y,vals,labI]=mu_util('ygrid',MAP_VAR_LIST.lats,MAP_VAR_LIST.longs,varargin{1},31,varargin{2:3});

  case 'box'

    [X,Y]=mu_util('box',31);

end

author: goGPS-Project


text:
function [gal] = ft32gal(ft3)
% Convert volume from cubic feet to US liquid gallons.
% Chad Greene 2012
gal = ft3*7.4805194805;

author: Sable

repo: mcbench-benchmarks


text:
function this = load(this, inputDataOrFile, varargin)
% loads (meta-)data from file(s), order defined by loopDimensions
%
%   Y = MrDataNd()
%   Y.load(varargin)
%
% This is a method of class MrDataNd.
%
% IN
%   inputDataOrFile     can be one of the following inputs
%                       1)  a Matrix: MrDataNd is created along with a
%                           dimInfo matching the dimensions of the data
%                           matrix
%                       2)  a figure/axes handle: MrDataNd tries to infer
%                           from the plot type and data, what Image shall
%                           be loaded from the specified handle
%                       3)  a file-name: MrDataNd is loaded from the
%                           specified file
%                       4)  cell(nFiles,1) of file names to be concatenated
%                       5)  a directory: All image files in the specified
%                           directory
%                       6)  a regular expression for all file names to be
%                           selected
%                           e.g. 'folder/fmri.*\.nii' for all nifti-files
%                           in a folder
%
%
%   varargin:   propertyName/value pairs, referring to
%               a) loading of files, e.g. 'updateProperties' or
%               'selectedVolumes'
%               b) 'select' struct to select a subset of data
%               c) 'dimInfo' object
%               d) property/value pairs for dimInfo
%               e) 'affineTrafo' object
%
%
%
% OUT
%   this        MrDataNd with updated .data and .dimInfo
%   affineTransformation
%               For certain file types, the affineTransformation is saved as a
%               header information. While ignored in MrDataNd, it might be
%               useful to return it for specific processing
%               See also MrImage MrAffineTransformation
%
% EXAMPLE
%   Y.load(gca)
%   Y.load(gcf);
%   Y.load(figure(121)); % to make distinction between fig handle and 1-element integer image
%
%   See also MrDataNd demo_save MrDataNd.read_matrix_from_workspace MrDataNd.read_data_from_graphics_handle

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2016-10-21
% Copyright (C) 2016 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% update/load path for data and dimInfo
% In load:
% --------
% Files are determined and loop over individual files is started.
%   In read_single_file:
% ----------------------
%   1:  Values are derived from the input matrix (nSamples) and/or file
%       (header info).
%       DimInfo is initiated here.
%   2:  If a _dimInfo.mat file exists, this is automatically loaded as well.
%       DimInfo properties are updated.
% End of read_single_file.
% ------------------------
% Single files are combined.
% 3:  If a dimInfo/affineTrafo object is an input argument,
%     dimInfo/affineTrafo properties are updated.
% 4:  If prop/val pairs are given,
%     dimInfo/affineTrafo properties are updated.

%% 0. Preliminaries
% process input parameters
if nargin < 2
    inputDataOrFile = this.get_filename();
end
defaults.select = [];
defaults.dimInfo = [];
defaults.affineTransformation = [];
[args, argsUnused] = tapas_uniqc_propval(varargin, defaults);
tapas_uniqc_strip_fields(args);

% Harvest propval for dimInfo constructor
[propValDimInfo, argsUnusedAfterDimInfo] = this.dimInfo.get_struct(argsUnused);

% input arguments for dimInfo constructor that are not properties
% themselves also need to be respected...
extraArgsDimInfo = this.dimInfo.get_additional_constructor_inputs();
[propValDimInfoExtra, argsUnusedAfterDimInfo] = tapas_uniqc_propval(argsUnusedAfterDimInfo, extraArgsDimInfo);

% create propValAffineTrafo
affineTrafo = MrAffineTransformation();
[propValAffineTransformation, loadInputArgs] = affineTrafo.get_struct(argsUnusedAfterDimInfo);
% check inputs
hasInputDimInfo = ~isempty(dimInfo);
hasInputAffineTransformation = ~isempty(affineTransformation);
hasPropValDimInfo = any(structfun(@(x) ~isempty(x), propValDimInfo));
hasPropValDimInfoExtra = any(structfun(@(x) ~isempty(x), propValDimInfoExtra));
hasPropValAffineTransformation = any(structfun(@(x) ~isempty(x), propValAffineTransformation));

hasSelect = ~isempty(select);
doLoad = 1;
%% 1. Determine files (for wildcards or folders)
isMatrix = isnumeric(inputDataOrFile) || islogical(inputDataOrFile);
isFigureOrAxesHandle = isa(inputDataOrFile, 'matlab.ui.Figure') || ...
    isa(inputDataOrFile, 'matlab.graphics.axis.Axes');
if isMatrix
    this.read_matrix_from_workspace(inputDataOrFile);
elseif isFigureOrAxesHandle
    this.read_data_from_graphics_handle(inputDataOrFile);
else % files or file pattern or directory
    isExplicitFileArray = iscell(inputDataOrFile) && ischar(inputDataOrFile{1});

    if isExplicitFileArray
        fileArray = inputDataOrFile;
    else
        fileArray = tapas_uniqc_get_filenames(inputDataOrFile);
    end

    % remove _dimInfo.mat from fileArray list
    [~, fileArray] = tapas_uniqc_find_info_file(fileArray, '_dimInfo.mat');
    %  get extra dimInfos from file names for select
    dimInfoExtra = MrDimInfo();
    dimInfoExtra.set_from_filenames(fileArray);
    % remove singleton dimensions
    dimInfoExtra.remove_dims();
    % now use select to only load subset of files
    % split select to within file (x, y, z, t) or between files
    [~, selectBetweenFiles, selectInFile] = dimInfoExtra.select(select);

    % number of files
    nFiles = numel(fileArray);

    %% 2. Load individual files into array of MrDataNd (including data of MrDimInfo)
    if nFiles == 1
        % only one file, read_single_files does everything that's necessary
        imgRead = read_single_file(this, fileArray{1}, loadInputArgs{:});
        % apply general select here (no splitting necessary)
        this.update_properties_from(imgRead.select(select));
    else
        % loop over nFiles and load each individually
        % initialize dataNdArray
        dataNdArray = cell(nFiles, 1);

        % actual constructor of possible sub-class of MrDataNd is
        % used here in loop, to reuse load as is in subclasses
        handleClassConstructor = str2func(class(this));

        for iFile = 1:nFiles
            % get filename
            fileName = fileArray{iFile};
            % get dimLabels from file name
            [dimLabels, dimValues] = tapas_uniqc_get_dim_labels_from_string(fileName);
            % check if dimLabels could be inferred from filename
            hasFoundDimLabelInFileName = ~isempty(dimLabels);
            if hasSelect
                % check if dimLabels and dimValues of this file are part of
                % selectDimInfo, to avoid unnecessary loading
                doLoad = all(cellfun(@ismember, num2cell(dimValues), selectBetweenFiles'));
            end
            if doLoad
                % load file into new file
                fprintf('Loading File %d/%d\n', iFile, nFiles);

                dataNdArray{iFile} = handleClassConstructor(fileName, ...
                    'select', selectInFile, loadInputArgs{:});
                % generate additional dimInfo
                if hasFoundDimLabelInFileName
                    % add units as samples
                    [units(1:numel(dimLabels))] = {'samples'};
                else
                    % generate generic dimLabels
                    dimLabels = {'file'};
                    dimValues = iFile;
                    units = 'sample';
                end
                % check if dimLabels already read
                hasDimLabel = any(ismember(dimLabels, dataNdArray{iFile}.dimInfo.dimLabels));
                if ~hasDimLabel
                    % add dimLabel and dim Value
                    dimsToAdd = dataNdArray{iFile}.dimInfo.nDims+1:dataNdArray{iFile}.dimInfo.nDims+numel(dimLabels);
                    dataNdArray{iFile}.dimInfo.add_dims(dimsToAdd, ...
                        'dimLabels', dimLabels, 'samplingPoints', dimValues,...
                        'units', units);
                end
            end
        end
        %% 3. Use combine to create one object
        % remove all empty cells from dataNdArray
        dataNdArray(cellfun(@isempty, dataNdArray)) = [];
        % use combine to create composite image
        imagesCombined = dataNdArray{1}.combine(dataNdArray);
        % add data to this
        this.update_properties_from(imagesCombined);
    end
end
% update dimInfo using input dimInfo
if hasInputDimInfo
    this.dimInfo.update_and_validate_properties_from(dimInfo);
end

% update dimInfo using prop/val dimInfo
if hasPropValDimInfo
    this.dimInfo.update_and_validate_properties_from(propValDimInfo);
end

if hasPropValDimInfoExtra
    this.dimInfo.set_dims(1:this.dimInfo.nDims, propValDimInfoExtra);
end

% update affineTransformation using input affineTransformation
if hasInputAffineTransformation
    this.affineTransformation.update_properties_from(affineTransformation);
end

% update affineTransformation using prop/val affineTransformation
if hasPropValAffineTransformation
    % affine matrix is a dependent property - cannot be changed via
    % update_properties_from
    if ~isempty(propValAffineTransformation.affineMatrix)
        this.affineTransformation.update_from_affine_matrix(propValAffineTransformation.affineMatrix)
    end
    this.affineTransformation.update_properties_from(propValAffineTransformation);
end
end

author: translationalneuromodeling

repo: tapas


text:
%This function makes a simple avi movie from a sequence of frames
%The user can control the display time of each frame. The movie is created
%in the same folder where this function is run.
%
%Usage:
%      Inputs:
%             name: root name of the framse
%         filetype: '.bmp'   ,  '.jpg'  , '.jpeg' , '.tif'  ,,,,
% number_of_frames: number of the frames
% display_time_of_frame: it is the display time of each frame
%
%           Output:
%  it creates a file named as Movie.avi in the same folder of this function
%
%            Usage: Try the simple provided example
%   Movie_from_frames('image','.bmp',4,10)
%This function is written by :
%                             Nassim Khaled
%                             American University of Beirut
%                             Research Assistant and Masters Student

function Movie_from_frames(name,filetype,number_of_frames,display_time_of_frame)
mov = avifile('MOVIE.avi');
count=0;
for i=1:number_of_frames
    name1=strcat(name,num2str(i),filetype);
    a=imread(name1);
    while count<display_time_of_frame
        count=count+1;
        imshow(a);
        F=getframe(gca);
        mov=addframe(mov,F);
    end
    count=0;
end
close all
mov=close(mov);


text:
function [expr morder] = sim_ex_Mullen_2011_Seizure(varargin)
% Simulation:  Epileptic Seizure
%
% Description:
%
%   13-variate VAR[6] system of stocastically-forced, damped coupled oscillators with time-varying (non-stationary) coupling.
%
%   This simulation creates a simulated "seizure" with time-varying coupling between clusters of sources which switches between 4 different stages.
%
%   The simulation is designed to be single-trial.
%
%   Recommended Settings:
%   Number of trials: 1
%   Sampling Rate:    100
%
%   The directed graph for this model can be viewed by executing the following command:
%
%   >>hlp_viewGraphicsResource('/sim/Mullen_2011_Seizure.jpg');
%
%
% Author Credits:
%
%   Tim Mullen, 2011
%
% References and Code:
%
%   Mullen, 2011. Unpublished
%
% ------------------------------------------------------------------------

% Cluster 1
f1   = .2;      % central frequency
tau1 = .2;      % damping time (larger-->highest variance)
% Cluster 2
f2   = .20;
tau2 = 3;
% Cluster 3
f3   = .10;
tau3 = 7;
% Cluster 4
f4   = .10;
tau4 = 6;

% start of seizure (1 minute)
S0 = 6000;
% set the approximate durations of each stage
S1_duration = 500;  % 5 sec at 100 Hz
S2_duration = 500;
S3_duration = 500;
S4_duration = 500;
% set seizure stage mid-points
S1_center = S0+S1_duration/2;                      % center of stage 1
S2_center = S1_center+S1_duration;                 % center of stage 2
S3_center = S2_center;                             % center of stage 3
S4_center = S3_center+S3_duration;                 % center of stage 4

Offset = 0;
SamplingRate = 1;

% specify the default system of equations
expr_def = { ...
          sprintf('x1(t) = {2*exp(-1/(%f+normpdfg(t,%f,%f,%f,100)))*cos(2*pi*%f/1)}*x1(t-1) + {-exp(-2/(%f+normpdfg(t,%f,%f,%f,100)))}*x1(t-2) + e1(t)',tau1,S1_duration/2,8,Offset+S1_center, f1, tau1,S1_duration,8,Offset+S1_center) ...                    % Ictal driver
          sprintf('x2(t) = %s + {normpdfg(t,%f,%f,%f,0.01)}*x3(t-2) + {normpdfg(t,%f,%f,%f,0.01)}*x4(t-2)  + {normpdfg(t,%f,%f,%f,1.3)}*x1(t-6)    + e2(t)',sim_dampedOscillator(f2,tau2,SamplingRate,2),  S1_duration/2,8,Offset+S1_center,   S1_duration/2,8,Offset+S1_center,  S1_duration/2,8,Offset+S1_center) ...                     % CLUSTER 2
          sprintf('x3(t) = %s +  {normpdfg(t,%f,%f,%f,0.7)}*x2(t-2) + {normpdfg(t,%f,%f,%f,0.1)}*x4(t-2)  + {normpdfg(t,%f,%f,%f,0.3)}*x5(t-3)    + e3(t)',sim_dampedOscillator(f2-1,tau2,SamplingRate,3),  S1_duration/2,8,Offset+S1_center,   S1_duration/2,8,Offset+S1_center,   (S2_duration+S3_duration)/2,10,Offset+S3_center) ...       % CLUSTER 2
          sprintf('x4(t) = %s +  {normpdfg(t,%f,%f,%f,0.7)}*x2(t-2) + {normpdfg(t,%f,%f,%f,0.1)}*x3(t-2)                                          + e4(t)',sim_dampedOscillator(f2+1,tau2,SamplingRate,4),  S1_duration/2,8,Offset+S1_center,  S1_duration/2,8,Offset+S1_center) ... % CLUSTER 2
          sprintf('x5(t) = %s +  {normpdfg(t,%f,%f,%f,0.001)}*x6(t-2) + {normpdfg(t,%f,%f,%f,0.5)}*x3(t-3)        + e5(t)'  ,sim_dampedOscillator(f3,tau3,SamplingRate,5), S3_duration/2,8,Offset+S3_center , S3_duration/2,8,Offset+S3_center) ...  % CLUSTER 3
          sprintf('x6(t) = %s +  {normpdfg(t,%f,%f,%f,0.001)}*x5(t-2)                                             + e6(t)'  ,sim_dampedOscillator(f3,tau3,SamplingRate,6), S3_duration/2,8,Offset+S3_center) ...                                  % CLUSTER 3
          sprintf('x7(t) = %s + {normpdfg(t,%f,%f,%f,1.3)}*x4(t-6) + {normpdfg(t,%f,%f,%f,0.01)}*x9(t-2) + {normpdfg(t,%f,%f,%f,0.01)}*x8(t-2) + {normpdfg(t,%f,%f,%f,0.01)}*x10(t-2)    + e7(t)' ,sim_dampedOscillator(f4,tau4,SamplingRate,7),     S4_duration/2,8,Offset+S4_center,     S4_duration/2,8,Offset+S4_center,     S4_duration/2,8,Offset+S4_center,     S4_duration/2,8,Offset+S4_center) ...  % CLUSTER 4
          sprintf('x8(t) = %s + {normpdfg(t,%f,%f,%f,0.6)}*x7(t-2)         + e8(t)' ,sim_dampedOscillator(f4,tau4,SamplingRate,8),     S4_duration/2,8,Offset+S4_center) ... % CLUSTER 4
          sprintf('x9(t) = %s + {normpdfg(t,%f,%f,%f,0.6)}*x7(t-2)         + e9(t)' ,sim_dampedOscillator(f4+1,tau4,SamplingRate,9),   S4_duration/2,8,Offset+S4_center) ... % CLUSTER 4
          sprintf('x10(t)= %s + {normpdfg(t,%f,%f,%f,0.6)}*x7(t-2)         + e10(t)',sim_dampedOscillator(f4-1,tau4,SamplingRate,10),  S4_duration/2,8,Offset+S4_center) ... % CLUSTER 4
          sprintf('x11(t)= 1.3*x11(t-1) + -0.8*x11(t-2)                             + e11(t)') ...    % Decoupled
          sprintf('x12(t)= 1.2*x12(t-1) + -0.4*x12(t-2)                             + e12(t)') ...    % Decoupled
          sprintf('x13(t)= 0.8*x13(t-1) + -0.4*x13(t-2) + -0.1*x13(t-4)             + e13(t)') ...    % Decoupled
          };

% set up argument definitions
arg_define(varargin, ...
    arg({'expr','DynamicalEquations'},expr_def,[],'System of equations'), ...
    arg({'morder','ModelOrder'},6,[1 Inf],'Model order. This is mandatory'), ...
    arg_subtoggle({'setDynamics','SetDynamics'},'off', ...
        {...
            arg_sub({'c1','Cluster1'},{},...
            { ...
                arg({'f0','Freq'},f1,[0 Inf],'Central frequency (Hz). This is in normalized Hz (i.e. f0=0.5=SamplingRate/2).'), ...
                arg({'tau','DampingTime'},tau1,[0 Inf],'Damping time') ...
            },'Cluster 1 Settings'), ...
            arg_sub({'c2','Cluster2'},{},...
            { ...
                arg({'f0','Freq'},f2,[0 Inf],'Central frequency (Hz). This is in normalized Hz (i.e. f0=0.5=SamplingRate/2).'), ...
                arg({'tau','DampingTime'},tau2,[0 Inf],'Damping time') ...
            },'Cluster 2 Settings'), ...
            arg_sub({'c3','Cluster3'},{},...
            { ...
                arg({'f0','Freq'},f3,[0 Inf],'Central frequency (Hz). This is in normalized Hz (i.e. f0=0.5=SamplingRate/2).'), ...
                arg({'tau','DampingTime'},tau3,[0 Inf],'Damping time') ...
            },'Cluster 3 Settings'), ...
            arg_sub({'c4','Cluster4'},{},...
            { ...
                arg({'f0','Freq'},f4,[0 Inf],'Central frequency (Hz). This is in normalized Hz (i.e. f0=0.5=SamplingRate/2).'), ...
                arg({'tau','DampingTime'},tau4,[0 Inf],'Damping time') ...
            },'Cluster 4 Settings'), ...
            arg_sub({'s1','Stage1'},{},...
            { ...
                arg({'duration','Duration'},S1_duration,[0 Inf],'Duration (samples)'), ...
                arg({'center','Center'},S1_center,[0 Inf],'Midpoint of stage (samples)'), ...
            },'Stage 1 Settings'), ...
            arg_sub({'s2','Stage2'},{},...
            { ...
                arg({'duration','Duration'},S2_duration,[0 Inf],'Duration (samples)'), ...
                arg({'center','Center'},S2_center,[0 Inf],'Midpoint of stage (samples)'), ...
            },'Stage 2 Settings'), ...
            arg_sub({'s3','Stage3'},{},...
            { ...
                arg({'duration','Duration'},S3_duration,[0 Inf],'Duration (samples)'), ...
                arg({'center','Center'},S3_center,[0 Inf],'Midpoint of stage (samples)'), ...
            },'Stage 3 Settings'), ...
            arg_sub({'s4','Stage4'},{},...
            { ...
                arg({'duration','Duration'},S4_duration,[0 Inf],'Duration (samples)'), ...
                arg({'center','Center'},S4_center,[0 Inf],'Midpoint of stage (samples)'), ...
            },'Stage 4 Settings'), ...
            arg({'szStart','SeizureStart'},S0,[0 Inf],'Start of seizure (samples)'), ...
        },'Override default system dynamics'));

if isempty(morder)
    error('SIFT:sim_examples:badParam','ModelOrder must be specified');
end

if g.setDynamics.arg_selection
    % re-evaluate expression with user-defined settings

    % set central frequency and damping times
    % Cluster 1
    f1   = g.setDynamics.c1.f0;      % central frequency
    tau1 = g.setDyamics.c1.tau;      % damping time (larger-->highest variance)
    % Cluster 2
    f2   = g.setDynamics.c2.f0;
    tau2 = g.setDyamics.c2.tau;
    % Cluster 3
    f3   = g.setDyamics.c3.f0;
    tau3 = g.setDyamics.c3.tau;
    % Cluster 4
    f4   = g.setDyamics.c4.f0;
    tau4 = g.setDyamics.c4.tau;

    % set start of seizure
    S0 = g.setDynamics.szStart;

    % set the approximate durations of each stage
    S1_duration = g.setDynamics.s1.duration;
    S2_duration = g.setDynamics.s2.duration;
    S3_duration = g.setDynamics.s3.duration;
    S4_duration = g.setDynamics.s4.duration;

    % set seizure stage mid-points
    if ~isempty(g.setDynamics.s1.center)
        S1_center = g.setDynamics.s1.center;
    else
        S1_center = S0+S1_duration/2;
    end
    if ~isempty(g.setDynamics.s2.center)
        S2_center = g.setDynamics.s2.center;
    else
        S2_center = S1_center+S1_duration;
    end
    if ~isempty(g.setDynamics.s3.center)
        S3_center = g.setDynamics.s3.center;
    else
         S3_center = S2_center;
    end
    if ~isempty(g.setDynamics.s4.center)
        S4_center = g.setDynamics.s4.center;
    else
        S4_center = S3_center+S3_duration;
    end

    % generate system of equations
    expr = { ...
          sprintf('x1(t) = {2*exp(-1/(%f+normpdfg(t,%f,%f,%f,100)))*cos(2*pi*%f/1)}*x1(t-1) + {-exp(-2/(%f+normpdfg(t,%f,%f,%f,100)))}*x1(t-2) + e1(t)',tau1,S1_duration/2,8,Offset+S1_center, f1, tau1,S1_duration,8,Offset+S1_center) ...                    % Ictal driver
          sprintf('x2(t) = %s + {normpdfg(t,%f,%f,%f,0.01)}*x3(t-2) + {normpdfg(t,%f,%f,%f,0.01)}*x4(t-2)  + {normpdfg(t,%f,%f,%f,1.3)}*x1(t-6)    + e2(t)',sim_dampedOscillator(f2,tau2,SamplingRate,2),  S1_duration/2,8,Offset+S1_center,   S1_duration/2,8,Offset+S1_center,  S1_duration/2,8,Offset+S1_center) ...                     % CLUSTER 2
          sprintf('x3(t) = %s +  {normpdfg(t,%f,%f,%f,0.7)}*x2(t-2) + {normpdfg(t,%f,%f,%f,0.1)}*x4(t-2)  + {normpdfg(t,%f,%f,%f,0.3)}*x5(t-3)    + e3(t)',sim_dampedOscillator(f2-1,tau2,SamplingRate,3),  S1_duration/2,8,Offset+S1_center,   S1_duration/2,8,Offset+S1_center,   (S2_duration+S3_duration)/2,10,Offset+S3_center) ...       % CLUSTER 2
          sprintf('x4(t) = %s +  {normpdfg(t,%f,%f,%f,0.7)}*x2(t-2) + {normpdfg(t,%f,%f,%f,0.1)}*x3(t-2)                                          + e4(t)',sim_dampedOscillator(f2+1,tau2,SamplingRate,4),  S1_duration/2,8,Offset+S1_center,  S1_duration/2,8,Offset+S1_center) ... % CLUSTER 2
          sprintf('x5(t) = %s +  {normpdfg(t,%f,%f,%f,0.001)}*x6(t-2) + {normpdfg(t,%f,%f,%f,0.5)}*x3(t-3)        + e5(t)'  ,sim_dampedOscillator(f3,tau3,SamplingRate,5), S3_duration/2,8,Offset+S3_center , S3_duration/2,8,Offset+S3_center) ...  % CLUSTER 3
          sprintf('x6(t) = %s +  {normpdfg(t,%f,%f,%f,0.001)}*x5(t-2)                                             + e6(t)'  ,sim_dampedOscillator(f3,tau3,SamplingRate,6), S3_duration/2,8,Offset+S3_center) ...                                  % CLUSTER 3
          sprintf('x7(t) = %s + {normpdfg(t,%f,%f,%f,1.3)}*x4(t-6) + {normpdfg(t,%f,%f,%f,0.01)}*x9(t-2) + {normpdfg(t,%f,%f,%f,0.01)}*x8(t-2) + {normpdfg(t,%f,%f,%f,0.01)}*x10(t-2)    + e7(t)' ,sim_dampedOscillator(f4,tau4,SamplingRate,7),     S4_duration/2,8,Offset+S4_center,     S4_duration/2,8,Offset+S4_center,     S4_duration/2,8,Offset+S4_center,     S4_duration/2,8,Offset+S4_center) ...  % CLUSTER 4
          sprintf('x8(t) = %s + {normpdfg(t,%f,%f,%f,0.6)}*x7(t-2)         + e8(t)' ,sim_dampedOscillator(f4,tau4,SamplingRate,8),     S4_duration/2,8,Offset+S4_center) ... % CLUSTER 4
          sprintf('x9(t) = %s + {normpdfg(t,%f,%f,%f,0.6)}*x7(t-2)         + e9(t)' ,sim_dampedOscillator(f4+1,tau4,SamplingRate,9),   S4_duration/2,8,Offset+S4_center) ... % CLUSTER 4
          sprintf('x10(t)= %s + {normpdfg(t,%f,%f,%f,0.6)}*x7(t-2)         + e10(t)',sim_dampedOscillator(f4-1,tau4,SamplingRate,10),  S4_duration/2,8,Offset+S4_center) ... % CLUSTER 4
          sprintf('x11(t)= 1.3*x11(t-1) + -0.8*x11(t-2)                             + e11(t)') ...    % Decoupled
          sprintf('x12(t)= 1.2*x12(t-1) + -0.4*x12(t-2)                             + e12(t)') ...    % Decoupled
          sprintf('x13(t)= 0.8*x13(t-1) + -0.4*x13(t-2) + -0.1*x13(t-4)             + e13(t)') ...    % Decoupled
          };
end

repo: REST


text:
classdef recovery_problems

    methods(Static)

        function problem = problem_small_1()
            m = 64;
            n = 121;
            k = 4;
            dict = spx.dict.simple.gaussian_dict(m, n);
            gen = spx.data.synthetic.SparseSignalGenerator(n, k);
            % create a sparse vector
            rep =  gen.biGaussian();
            signal = dict*rep;
            problem.dictionary = dict;
            problem.representation_vector = rep;
            problem.sparsity_level = k;
            problem.signal_vector = signal;
        end

        function problem = problem_large_1()
            m = 1000;
            n = 2000;
            k = 100;
            dict = spx.dict.simple.gaussian_dict(m, n);
            gen = spx.data.synthetic.SparseSignalGenerator(n, k);
            % create a sparse vector
            rep =  gen.biGaussian();
            signal = dict*rep;
            problem.dictionary = dict;
            problem.representation_vector = rep;
            problem.sparsity_level = k;
            problem.signal_vector = signal;
        end

        function problem = problem_barbara_blocks()
            env = spx_get_env();
            images_dir = env.local_settings.standard_test_images_dir;
            image_path = fullfile(images_dir, 'barbara.png');
            image = spx.data.synthetic.recovery_problems.read_image(image_path);
            blkSize = 8;
            patches = im2col(image, [blkSize, blkSize], 'distinct');
            problem.signals = patches;
            problem.image = image;
            problem.blkSize = blkSize;
            problem.sparsity_level = 4;
            % Let's prepare a dictionary for the image
            N = 64;
            D = 121;
            problem.dictionary = spx.dict.simple.overcomplete2DDCT(N, D);
        end

        function problem = problem_test_image_blocks(image_name, block_type)
            if nargin < 2
                block_type = 'distinct';
            end
            env = spx_get_env();
            rootdir = env.local_settings.standard_test_images_dir;
            switch image_name
                case 'barbara'
                    fname = 'barbara.png';
                case 'cameraman'
                    fname = 'cameraman.tif';
                case 'house'
                    fname = 'house.tif';
                case 'jetplane'
                    fname = 'jetplane.tif';
                case 'lake'
                    fname = 'lake.tif';
                case 'lena'
                    fname = 'lena_gray_512.tif';
                case 'livingroom'
                    fname = 'livingroom.tif';
                case 'mandril'
                    fname = 'mandril_gray.tif';
                case 'peppers'
                    fname = 'peppers_gray.tif';
                case 'pirate'
                    fname = 'pirate.tif';
                case 'walkbridge'
                    fname = 'walkbridge.tif';
                case 'blonde'
                    fname = 'woman_blonde.tif';
                case 'darkhair'
                    fname = 'woman_darkhair.tif';
                otherwise
                    error('Unsupported test image');
            end
            filepath = fullfile(rootdir, fname);
            image = spx.data.synthetic.recovery_problems.read_image(filepath);
            blkSize = 8;
            patches = im2col(image, [blkSize, blkSize], block_type);
            problem.signals = patches;
            problem.image = image;
            problem.blkSize = blkSize;
        end

        function image = read_image(filepath)
            [image,~]=imread(filepath);
            image = im2double(image);
            % If the image is RGB, let us convert it to gray
            sz = size(image);
            if (length(sz) > 2)
                if sz(3) == 2
                    image = image(:, :, 1);
                elseif sz(3) == 3
                    image = rgb2gray(image);
                end
            end
            % Let us bring image to 0-255 range
            if max(image(:)) < 2
                image = image * 255;
            end
        end

    end

end

author: indigits

repo: sparse-plex


text:
function rs = realifft(data, N)

% inverse fft for fourier coefficents from real valued  original data
% needs the length of the original time series from which the fft was computed
% see also : realfft

data = data(:);
if mod(N,2) 	% odd length
	data = [data ; conj(data(end:-1:2))];
else			% even length
	data = [data ; conj(data(end-1:-1:2))];
end

rs = real(ifft(data)); % remove small imaginary components introduced by rounding errors

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function Write_Kinematics(FileWritingSetup)
% Write_Kinematics(FileWritingSetup)
%
% FUNCTION:
%
%   This function takes the symbolic expressions for system energy,
%   position, and velocity, and uses them to write a function file.
%
%
% INPUTS: FileWritingSetup.
%
%   States = [2N x 2] cell array where the first column is the name of
%       the states, and the second column is comments about that states.
%       The index of a state in this array will match the index in the
%       dynamics function to be written. This function assumes a second
%       order system that is being written in first order form, so if
%       States{i,1} = z, then States{i+N,1} = dz; Furthermore, the output
%       (dStates) will have the form: dStates(i,:) = dz,
%       and dStates(i+N,:) = ddz.
%
%   Parameters = [K x 2] cell array with all of the input parameters. The
%       first column is for names, the second for comments
%
%   Kinematics = A struct with fields for the symbolic expressions for
%       position and velocity.
%
%   Directory = the name of the directory to write the files in.
%
%
% OUTPUTS:
%   This function will write the function kinematics()
%
%
% NOTE:
%
%   1)  It is assumed that the symbols in all of the inputs are self
%       consistent; that is - if Contacts{1,1} = H1, then there must be a field
%       in Dyn.(FileNames{i,1}) called H1.
%
%
% Written by Matthew Kelly
% December 6, 2013
% Cornell University
%
% See also DERIVE_EOM

States = FileWritingSetup.States;
Kinematics = FileWritingSetup.Kinematics;
Parameters = FileWritingSetup.Parameters;
Directory = FileWritingSetup.Directory;
Energy = FileWritingSetup.Energy;

r0 = Kinematics.r0;
r1 = Kinematics.r1;
r2 = Kinematics.r2;
dr0 = Kinematics.dr0;
dr1 = Kinematics.dr1;
dr2 = Kinematics.dr2;

CurrDir = cd;
if ~exist(Directory,'dir')
    mkdir(Directory);
end
cd(Directory);

fid = fopen('kinematics.m','w');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                         Function Header                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

fprintf(fid,['function [Position, Velocity, Energy] = ' ...
    'kinematics(States, Parameters)\n']);
fprintf(fid,['%% function [Position, Velocity, Energy] = ' ...
    'kinematics(States, Parameters)\n']);
fprintf(fid,'%%\n%% Computer Generated File -- DO NOT EDIT \n');
fprintf(fid,['%%\n%% This function was created by the function '...
    'Write_Kinematics()\n']);
fprintf(fid,['%% ' datestr(now) '\n%%\n']);
fprintf(fid,'%%  RETURNS: \n');
fprintf(fid,'%%   Energy = a struct with fields: \n');
fprintf(fid,'%%       Potential \n');
fprintf(fid,'%%       Kinetic \n');
fprintf(fid,'%%       Total \n');
fprintf(fid,'%% \n');
fprintf(fid,'%%   Position = a struct with fields: \n');
fprintf(fid,'%%       footOne = [2xN] matrix with the position vectors for Foot One \n');
fprintf(fid,'%%       footTwo = [2xN] matrix with the position vectors for Foot Two \n');
fprintf(fid,'%%       hip = [2xN] matrix with the position vectors for the hip \n');
fprintf(fid,'%%       CoM = [2xN] matrix with the position vectors for the center of mass\n');
fprintf(fid,'%% \n');
fprintf(fid,'%%   Velocity = a struct with fields: \n');
fprintf(fid,'%%       footOne = [2xN] matrix with the velocity vectors for Foot One \n');
fprintf(fid,'%%       footTwo = [2xN] matrix with the velocity vectors for Foot Two \n');
fprintf(fid,'%%       hip = [2xN] matrix with the velocity vectors for the hip \n');
fprintf(fid,'%%       CoM = [2xN] matrix with the velocity vectors for the center of mass \n');
fprintf(fid,'%%\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% Matthew Kelly \n%% Cornell University \n%% \n\n');
fprintf(fid,'%% \n');
fprintf(fid,'%% See also DERIVE_EOM \n');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                           Parse Inputs                              %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

for idxState=1:size(States,1)
    fprintf(fid,[States{idxState,1} ' = States(' num2str(idxState) ',:); %% ' ...
        States{idxState,2} '\n']);
end %idxState
fprintf(fid,'\n');

for idxParam=1:size(Parameters,1)
    fprintf(fid,[Parameters{idxParam,1} ' = Parameters.' ...
        Parameters{idxParam,1} '; %% ' ...
        Parameters{idxParam,2} '\n']);
end %idxState
fprintf(fid,'\n');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                          Write Outputs                              %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

fprintf(fid,'N = size(States,2);\n');

fprintf(fid,'\n');
fprintf(fid,'Position.footOne = zeros(2,N);\n');
fprintf(fid,'Position.footTwo = zeros(2,N);\n');
fprintf(fid,'Position.hip = zeros(2,N);\n');
fprintf(fid,['Position.footOne(1,:) = ' vectorize(r0(1)) ';\n']);
fprintf(fid,['Position.footOne(2,:) = ' vectorize(r0(2)) ';\n']);
fprintf(fid,['Position.footTwo(1,:) = ' vectorize(r2(1)) ';\n']);
fprintf(fid,['Position.footTwo(2,:) = ' vectorize(r2(2)) ';\n']);
fprintf(fid,['Position.hip(1,:) = ' vectorize(r1(1)) ';\n']);
fprintf(fid,['Position.hip(2,:) = ' vectorize(r1(2)) ';\n']);
fprintf(fid,'Position.CoM = (m1*Position.footOne + m2*Position.footTwo + M*Position.hip)/(m1+m2+M);\n');
fprintf(fid,'\n');

fprintf(fid,'if nargout > 1');
fprintf(fid,'\n');
fprintf(fid,'    Velocity.footOne = zeros(2,N);\n');
fprintf(fid,'    Velocity.footTwo = zeros(2,N);\n');
fprintf(fid,'    Velocity.hip = zeros(2,N);\n');
fprintf(fid,['    Velocity.footOne(1,:) = ' vectorize(dr0(1)) ';\n']);
fprintf(fid,['    Velocity.footOne(2,:) = ' vectorize(dr0(2)) ';\n']);
fprintf(fid,['    Velocity.footTwo(1,:) = ' vectorize(dr2(1)) ';\n']);
fprintf(fid,['    Velocity.footTwo(2,:) = ' vectorize(dr2(2)) ';\n']);
fprintf(fid,['    Velocity.hip(1,:) = ' vectorize(dr1(1)) ';\n']);
fprintf(fid,['    Velocity.hip(2,:) = ' vectorize(dr1(2)) ';\n']);
fprintf(fid,'    Velocity.CoM = (m1*Velocity.footOne + m2*Velocity.footTwo + M*Velocity.hip)/(m1+m2+M);\n');
fprintf(fid,'\n');
fprintf(fid,'end\n');

fprintf(fid,'if nargout > 2');
fprintf(fid,'\n');
fprintf(fid,['    Energy.Potential.m1 = ' ...
    vectorize(Energy.Potential.m1) ';\n']);
fprintf(fid,['    Energy.Potential.m2 = ' ...
    vectorize(    Energy.Potential.m2) ';\n']);
fprintf(fid,['    Energy.Potential.M = ' ...
    vectorize(Energy.Potential.M) ';\n']);
fprintf(fid,'    Energy.Potential.Total = Energy.Potential.m1 + Energy.Potential.m2 + Energy.Potential.M;\n');
fprintf(fid,'\n');
fprintf(fid,['    Energy.Kinetic.m1 = ' ...
    vectorize(Energy.Kinetic.m1) ';\n']);
fprintf(fid,['    Energy.Kinetic.m2 = ' ...
    vectorize(Energy.Kinetic.m2) ';\n']);
fprintf(fid,['    Energy.Kinetic.M = ' ...
    vectorize(Energy.Kinetic.M) ';\n']);
fprintf(fid,'    Energy.Kinetic.Total = Energy.Kinetic.m1 + Energy.Kinetic.m2 + Energy.Kinetic.M;\n');
fprintf(fid,'\n');
fprintf(fid,'    Energy.Total = Energy.Kinetic.Total + Energy.Potential.Total;\n');
fprintf(fid,'\n');
fprintf(fid,'end\n');
fprintf(fid,'end\n');
fclose(fid);

cd(CurrDir);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  ~~~~~~~~~~~~~~~~~~~ Sub Functions  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function characterExpression = vectorize(symbolicExpression)

        characterExpression = char(symbolicExpression);
        characterExpression = strrep(characterExpression,'*','.*');
        characterExpression = strrep(characterExpression,'/','./');
        characterExpression = strrep(characterExpression,'^','.^');
end

author: MatthewPeterKelly

repo: dscTutorials


text:
classdef (HandleCompatible) DisplayNonScalarObjectAsTable < matlab.mixin.CustomDisplay
    % DisplayNonScalarObjectAsTable -
    % ---------------------------------------------------------------------
    % Abstract: This mixin class overrides displayNonScalarObject() to
    % display an array of objects as a table.
    %

    %   Copyright 2018-2019 The MathWorks, Inc.
    %
    % Auth/Revision:
    %   MathWorks Consulting
    %   $Author: rjackey $
    %   $Revision: 348 $  $Date: 2018-03-02 15:51:54 -0500 (Fri, 02 Mar 2018) $
    % ---------------------------------------------------------------------

    %% Public Methods
    methods (Sealed)

        function t = toDisplayTable(obj)
            % Convert the object to a table

            % Prepare indices as row names
            numObj = numel(obj);
            if isrow(obj)
                indices = string(1:numObj)';
            else
                varargout = cell(1,ndims(obj));
                [varargout{:}] = ind2sub(size(obj),1:numObj);
                indices = string(vertcat(varargout{:})');
                indices = join(indices,',',2);
            end
            rowNames = "(" + indices + ")";

            % Check for deleted handle objects
            if isa(obj,'handle')
                isDeleted = ~isvalid(obj);
                obj(isDeleted) = [];
                rowNames(isDeleted) = [];
            end

            % Gather property info
            props = string( properties(obj) )';

            % Preallocate a cell to store values
            numProps = numel(props);
            numObj = numel(obj);
            values = cell(numObj,numProps);

            % Populate cell with property values
            for pIdx = 1:numProps
                thisProp = props(pIdx);
                values(:,pIdx) = {obj.(thisProp)}';
            end %for

            % Convert to table
            t = cell2table(values,'VariableNames',props,'RowNames',rowNames);

        end %function

    end %methods

    %% Protected Methods
    methods (Sealed, Access=protected)

        function displayNonScalarObject(obj)

            % Format text to display
            className = matlab.mixin.CustomDisplay.getClassNameForHeader(obj);
            dimStr = matlab.mixin.CustomDisplay.convertDimensionsToString(obj);

            % Display the header
            if isa(obj,'matlab.mixin.Heterogeneous')
                fprintf('  %s Heterogeneous %s with common properties:\n\n',dimStr,className);
            else
                fprintf('  %s %s with properties:\n\n',dimStr,className);
            end

            % Show the group list in a table
            disp( obj.toDisplayTable() );

            if isa(obj,'handle') && any(~isvalid(obj))
               fprintf('  Object array contains deleted handles.\n');
            end

        end %function

    end %methods

end % classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [cum_ret, cumprod_ret, daily_ret, ra_ret, run_time]...
    = manager(strategy_name, dataset_name, varargins, opts)
% Manager:
% This program universally calls different OLPS strategies and datasets.
%
% function [cum_ret, cumprod_ret, daily_ret, ra_ret, run_time]...
%          = manager(strategy_name, dataset_name, varargins, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% ra_ret:  result analysis, output risk-adjusted return and other
% run_time: computational time for a algorithm on a dataset
%
% strategy_name: the string of a strategy.
% dataset_name: the string of a dataset.
% varargins: parameter array, depends on the algorithms
% opts: a variable to control the behavior of back testing. Described as
% follows.
%
% Option parameters:
%     opts.quiet_mode = 0/1;    display debug info[Y/N];
%     opts.display_interval = 500; display info time interval
%     opts.log_record = 0/1;    record the .log file[N/Y]
%     opts.mat_record = 0/1;    record the .mat file[N/Y]
%     opts.analyze_mode = 0/1;  analyze the algorithm? [N/Y]
%
% Example: [cum_ret, cumprod_ret, daily_ret, ra_ret, run_time] ...
%          = manager('market', 'nyse_o', {0}, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors: Doyen Sahoo
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Loading the Dataset stringed dataset_name
if (~opts.quiet_mode)
    fprintf(1, 'Running strategy %s on dataset %s\nLoading dataset %s.\n', strategy_name, dataset_name, dataset_name);
end
load(sprintf('../Data/%s.mat', dataset_name));
if (~opts.quiet_mode)
    fprintf(1, 'Finish loading dataset %s\nThe size of the dataset is %dx%d.\n', dataset_name, size(data, 1), size(data, 2));
end

%% Logging stage, format the log file name
if (opts.log_mode)
    dt = datestr(now, 'yyyy-mmdd-HH-MM-SS');
    file_name = ['../Log/Log/' strategy_name '-' dataset_name '-' dt '.txt'];
    fid=fopen(file_name, 'wt');

    fprintf(fid, 'Running strategy %s on dataset %s.\n', strategy_name, dataset_name);
else
    fid = 1;  % use 1 to ensure no leak
end

start_time = datestr(now, 'yyyy-mmdd-HH-MM-SS-FFF');

if (opts.log_mode)
    fprintf(fid, 'Start Time: %s.\n', start_time);
end
if (~opts.quiet_mode)
    fprintf(1, 'Start Time: %s.\n', start_time);
end

%%%%%%%%%%%%%%Calling Strategy Stage%%%%%%%%%%%%%%%%%%%
fprintf(1, '----Begin %s on %s-----\n', strategy_name, dataset_name);
start_watch = tic;

strategy_fun = [strategy_name '(fid, data, varargins, opts)'];
[cum_ret, cumprod_ret, daily_ret, daily_portfolio] = eval(strategy_fun);

tElapse = toc(start_watch);
fprintf(1, '----End %s on %s-----\n', strategy_name, dataset_name);

isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
if (isOctave)
    fflush(stdout); % force output in Octave
end

%%%%%%%%%%%%%%Logging Stage%%%%%%%%%%%%%%%%%%%%%%
stop_time = datestr(now, 'yyyy-mmdd-HH-MM-SS-FFF');

if (opts.log_mode)
    fprintf(fid, 'Stop Time: %s.\n', stop_time);
    fprintf(fid, 'Elapse time(s): %f.\n', tElapse);
end

if (~opts.quiet_mode)
    fprintf(1, 'Stop Time: %s.\n', stop_time);
    fprintf(1, 'Elapse time(s): %f.\n', tElapse);
end

%%%%%%%%%%%%%%Analysis Stage%%%%%%%%%%%%%%%%%%%%%%
% Analyze the results, output to the log file and screen
ra_ret = [];
run_time = tElapse;
if (opts.analyze_mode)
    [ra_ret] = ra_result_analyze(fid, data, cum_ret, cumprod_ret, daily_ret, opts);
end

%% Save the result variables to a .mat file
if (opts.mat_mode)
    mat_dt = datestr(now, 'yyyy-mmdd-HH-MM');
    mat_name = ['../Log/Mat/' strategy_name '-' dataset_name '-' mat_dt '.mat'];
    save(mat_name,  'cum_ret', 'cumprod_ret', 'daily_ret', 'ra_ret', 'run_time', 'daily_portfolio');
end

% Clear the log file
if (opts.log_mode)
    fclose(fid);
end

%% Email notivation, currently disabled.
% %% Send the results to email: libin.yz@gmail.com
% if (opts.log_mode) && (opts.mat_mode)
%     msg_rcv = ['libin.yz@gmail.com'];
%     msg_head = ['Results from MATLAB:' strategy_name '-' dataset_name '-' dt ];
%     msg_body = ['cum_ret' num2str(cum_ret) '    ' 'run_time' num2str(run_time)];
%     msg_attach =[file_name];
% %    sendmail(msg_rcv, msg_head, msg_body,{msg_attach})
% end

end
%%%%%%%%%%%%%%End%%%%%%%%%%%%%%%%%%%%%%

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function PTKRemoveOrphanedCoronalSlices(mask_1, mask_2)
    % PTKRemoveOrphanedCoronalSlices.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    for coronal_index = 1 : mask_1.ImageSize(1);
        slice_1 = mask_1.GetSlice(coronal_index, PTKImageOrientation.Coronal);
        slice_2 = mask_2.GetSlice(coronal_index, PTKImageOrientation.Coronal);
        if any(slice_1(:)) && ~any(slice_2(:))
            slice_1(:) = 0;
            mask_1.ReplaceImageSlice(slice_1, coronal_index, PTKImageOrientation.Coronal);
        end
        if any(slice_2(:)) && ~any(slice_1(:))
            slice_2(:) = 0;
            mask_2.ReplaceImageSlice(slice_2, coronal_index, PTKImageOrientation.Coronal);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit