text:
function n = navier_q1(xy,ev,flowsol)
%navier_q1  Q1 convection matrix
%   N = navier_q1(xy,ev,flowsol);
%   input
%          xy         Q2 nodal coordinate vector
%          ev         element mapping matrix
%          flowsol    Q1-Q1 or Q1_P0 flow solution
%   output
%          N          Q1 scalar convection matrix
%
%   Natural boundary conditions apply.
%   IFISS function: DJS; 11 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
nngpt=4;
x=xy(:,1); y=xy(:,2);
nvtx=length(x);   nel=length(ev(:,1));
usol=flowsol(1:nvtx); vsol=flowsol(nvtx+1:2*nvtx);
fprintf('setting up Q1 convection matrix...  ')
%
% initialise global matrices
      n = sparse(nvtx,nvtx);
%
% Gauss point integration rules
      if (nngpt==4)        % 2x2 Gauss points
      gpt=1.0e0/sqrt(3.0e0);
      s(1) = -gpt; t(1) = -gpt; wt(1)=1;
      s(2) =  gpt; t(2) = -gpt; wt(2)=1;
      s(3) =  gpt; t(3) =  gpt; wt(3)=1;
      s(4) = -gpt; t(4) =  gpt; wt(4)=1;
      elseif (nngpt==1)   % 1x1 Gauss point
      s(1) =    0; t(1) =    0; wt(1)=4;
      else
	  error('Check Gauss point integration specification')
      end
%
% inner loop over elements
      for ivtx = 1:4
      xl_v(:,ivtx) = x(ev(:,ivtx));
      yl_v(:,ivtx) = y(ev(:,ivtx));
      xsl(:,ivtx) = usol(ev(:,ivtx));
	  ysl(:,ivtx) = vsol(ev(:,ivtx));
	  end
      ne = zeros(nel,4,4);
%
% loop over Gauss points
         for igpt = 1:nngpt
         sigpt=s(igpt);
         tigpt=t(igpt);
         wght=wt(igpt);
%  evaluate derivatives etc
         [jac,invjac,phi,dphidx,dphidy] = deriv(sigpt,tigpt,xl_v,yl_v);
         u_x = zeros(nel,1); u_y=zeros(nel,1);
            for k=1:4
		    u_x(:) = u_x(:) + xsl(:,k) .* phi(:,k);
		    u_y(:) = u_y(:) + ysl(:,k) .* phi(:,k);
		    end
		 for j = 1:4
            for i = 1:4
				ne(:,i,j)  = ne(:,i,j)  + wght*u_x(:).*phi(:,i).*dphidx(:,j);
                ne(:,i,j)  = ne(:,i,j)  + wght*u_y(:).*phi(:,i).*dphidy(:,j);
             end
	    end
%
% end of Gauss point loop
         end
%
%%  element assembly into global matrix
      for krow=1:4
	  nrow=ev(:,krow);
          for kcol=1:4
		  ncol=ev(:,kcol);
          n = n + sparse(nrow,ncol,ne(:,krow,kcol),nvtx,nvtx);
	      end
       end
%
fprintf('done.\n')
return


text:
function zD = cvx_newcnstr( zb, iseq, need_dual )

global cvx___
try
    pstr = cvx___.problems(end);
catch
    cvx_throw( 'No CVX model is present.' );
end
if isa( zb, 'cvx' ),
    zb = cvx_basis( zb );
end
[nN,mN] = size( zb );
isr = isreal( zb );
if ~isr,
    mN = 2 * mN;
    zb = reshape( [ real(zb) ; imag(zb) ], nN, mN );
end
if ~iseq,
    ndim = cvx_pushnneg( mN );
    zb(ndim(end),1) = 0;
    zb = zb + sparse( ndim, 1:mN, -1 );
end
cvx___.equalities{end+1} = zb;
cvx___.inequality(end+1) = ~iseq;
mO = cvx___.n_equality;
cvx___.n_equality = mO + mN;
cp = pstr.checkpoint(4);
if cp > 1,
    x = find( any( zb, 2 ), 2, 'first' );
    if any( x > 1 & x <= cp ),
        cvx___.problems(end).checkpoint(4) = x(1+(x(1)==1)) - 1;
    end
end
if nargout
    if nargin == 3 && need_dual
        zD = sparse( mO + 1 : mO + mN, 1 : mN, 1 );
        if ~isr, zD = zD(:,1:2:end) + 1j * zD(:,2:2:end); end
    else
        zD = [];
    end
end

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

author: yu-jiang

repo: radpbook

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


text:
%% Axis Distribution Function
%
%%

author: mtex-toolbox

repo: mtex

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


text:
classdef ThrottleWithNoThrustModelValidator < AbstractLaunchVehicleDataValidator
    %ThrottleWithNoThrustModelValidator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lvdData LvdData
    end

    methods
        function obj = ThrottleWithNoThrustModelValidator(lvdData)
            obj.lvdData = lvdData;
        end

        function [errors, warnings] = validate(obj)
            errors = LaunchVehicleDataValidationError.empty(0,1);
            warnings = LaunchVehicleDataValidationWarning.empty(0,1);

            warnEvtNums = [];
            evts = obj.lvdData.script.evts;
            for(i=1:length(evts))
                evt = evts(i);
                evtPropagator = evt.propagatorObj;
                propagatorEnum = evtPropagator.propagatorEnum;

                if(propagatorEnum == PropagatorEnum.ForceModel)
                    if(not(any(evtPropagator.forceModels == ForceModelsEnum.Thrust)))
                        stateLogEntries = obj.lvdData.stateLog.getAllStateLogEntriesForEvent(evt);
                        times = [stateLogEntries.time];
                        finalTime = max(times);

                        for(j=1:length(stateLogEntries))
                            stateLogEntry = stateLogEntries(j);
                            if(stateLogEntry.throttle > 0 && stateLogEntry.time < finalTime) %less than sign so we don't include the final time, which may have actions to enable the throttle and would be erronously picked up here
                                warnEvtNums(end+1) = evt.getEventNum(); %#ok<AGROW>
                                break;
                            end
                        end
                    end
                elseif(propagatorEnum == PropagatorEnum.TwoBody)
                    stateLogEntries = obj.lvdData.stateLog.getAllStateLogEntriesForEvent(evt);
                    times = [stateLogEntries.time];
                    finalTime = max(times);

                    for(j=1:length(stateLogEntries))
                        stateLogEntry = stateLogEntries(j);
                        if(stateLogEntry.throttle > 0 && stateLogEntry.time < finalTime) %less than sign so we don't include the final time, which may have actions to enable the throttle and would be erronously picked up here
                            warnEvtNums(end+1) = evt.getEventNum(); %#ok<AGROW>
                            break;
                        end
                    end
                end
            end

            if(not(isempty(warnEvtNums)))
                eventStr = makeEventsStr(unique(warnEvtNums));
                str = sprintf('Throttle is greater than 0%% on Events with Thrust model disabled. (Events: %s)\n%s', eventStr);
                warnings(end+1) = LaunchVehicleDataValidationWarning(str);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function res = gt(a,b)
%GT           Implements  a > b  elementwise for intervals a and b
%
%  if true,  a  is definitely greater than  b
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if ~isa(a,'intval')
    a = intval(a);
  end
  if ~isa(b,'intval')
    b = intval(b);
  end

  if a.complex | b.complex
    res = real(inf(a)) > real(sup(b)) & imag(inf(a)) > imag(sup(b)) ;
  else
    res = inf(a) > sup(b) ;
  end

repo: OpenMAS


text:
function fvsolver(N, method)
% N = number of grid points

% grid spacing on interval (-1,1)
dx = 2./N;

% location of cell centers
x = linspace(-1+0.5*dx, 1-0.5*dx, N);

% advection speed
u = 1;

% cfl condition + safety margin
dt = 0.2*dx;

% final time
FinalTime = 2;
Nsteps = ceil(FinalTime/dt);
dt = FinalTime/Nsteps;

% initial condition
q = fvexact(x); % q = 1*(abs(x)<=0.25);

% time step
for n=1:Nsteps

  qp1  = [q(2:N),q(1)];    % q_{m+1}
  qm1 = [q(N),q(1:N-1)]; % q_{m-1}

  % compute slopes s_m based on limiter choice

  if(strcmp(method,'minmod'))
      s = minmod((q-qm1)/dx, (qp1-q)/dx);
  elseif(strcmp(method,'MC'))
      s = minmod( (qp1-qm1)/(2*dx), minmod(2*(q-qm1)/dx, 2*(qp1-q)/dx));
  elseif(strcmp(method,'Fromm'))
      s = (qp1-qm1)/(2*dx);
  elseif(strcmp(method,'BeamWarming'))
      s = (q-qm1)/(dx); % upwind
  elseif(strcmp(method,'LaxWendroff'))
      s = (qp1-q)/(dx); % downwind
  else
      s = zeros(size(q));
  end

  sm1 = [s(N),s(1:N-1)];                           % s_{m-1}

  % update cell averages
  q = q - (u*dt/dx)*(q-qm1) - 0.5*(u*dt/dx)*(dx-u*dt)*(s-sm1);

  if(mod(n,40)==0)
    plot(x,q, '.'); pause(0.02);
  end
end

plot(x, q, 'b.');
% have gone one exact period
hold on; plot(x, fvexact(x), 'r-'); hold off;
xlabel('x');
legend(sprintf('%s N=%d T=%f', method, N, FinalTime), 'Exact solution')
axis([-1 1 -.5 1.5])

repo: Aero-matlab


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

J = sum(((X * Theta' - Y) .^ 2)(R == 1)) / 2;
J = J + lambda / 2 * (sum(sum(Theta .^ 2)) + sum(sum(X .^ 2)));

X_grad = ((X * Theta' - Y) .* R) * Theta;
Theta_grad = ((X * Theta' - Y) .* R)' * X;

X_grad = X_grad + (lambda * X);
Theta_grad = Theta_grad + (lambda * Theta);

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

grad = [X_grad(:); Theta_grad(:)];

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function varargout = chebpolyval3(f, varargin)
%CHEBPOLYVAL3   Values of a CHEBFUN3 object F on a tensor product grid.
%   X = CHEBPOLYVAL3(F, M, N, P) returns an M x N x P tensor of values of a
%   CHEBFUN3 object F on a tensor product grid.
%
%   X = CHEBPOLYVAL3(F) does the same but sets M, N, and P equal to the
%   length of F.
%
%   [CORE, C, R, T] = CHEBPOLYVAL3(F) returns the low rank representation
%   of the values of F on a tensor product grid, i.e.,
%   X = CORE x_1 C x_2 R x_3 T.
%
%   Example: For an order-3 discrete tensor R, we should have
%                                       R \approx chebpolyval3(chebfun3(R)).

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

% Empty check.
if ( isempty(f) )
    varargout = { [] };
    return
end

if ( nargin == 1 )
    % Get degrees:
    [m, n, p] = length(f);
elseif ( nargin ~= 4 )
    error('CHEBFUN:CHEBFUN3:chebpolyval3:inputs', 'Dimension not specified.');
else
    m = varargin{1};
    n = varargin{2};
    p = varargin{3};
end

% Get low rank representation of f:
[core, cols, rows, tubes] = tucker(f);

tech = chebfunpref().tech();
colVals = tech.coeffs2vals(chebcoeffs(cols, m));
rowVals = tech.coeffs2vals(chebcoeffs(rows, n));
tubeVals = tech.coeffs2vals(chebcoeffs(tubes, p));

% Evaluate:
if ( nargout <= 1 )
    varargout = {chebfun3.txm(chebfun3.txm(chebfun3.txm(core, colVals, 1), ...
        rowVals, 2), tubeVals, 3)};
else
    varargout = {core, colVals, rowVals, tubeVals};
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function value = f2 ( x )
  value = cos ( 3 * pi * x );
  return
end

author: johannesgerer

repo: jburkardt-m


text:
function[DepMat RgbMat] = Kinect()
%---------------------------------------------------------------
% This function retrive Dep and Rgb color matries from mex code
% And return as DepMat RgbMat, respectly, after reshape them.
%---------------------------------------------------------------
[a b] = getimagedata(2);
DepMat = reshape(a,640,480)';
temp = (reshape(b(1:3:end),640,480))';
RgbMat(:,:,1) = temp(:,end:-1:1);
temp = (reshape(b(2:3:end),640,480))';
RgbMat(:,:,2) = temp(:,end:-1:1);
temp = (reshape(b(3:3:end),640,480))';
RgbMat(:,:,3) = temp(:,end:-1:1);
end

author: Sable

repo: mcbench-benchmarks


text:
%-----------------------------------------------------------------------
% Job saved on 02-Jul-2014 11:27:33 by cfg_util (rev $Rev$)
% spm SPM - SPM12b (5672)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.spm.spatial.smooth.data = '<UNDEFINED>';
matlabbatch{1}.spm.spatial.smooth.fwhm = [8 8 8];
matlabbatch{1}.spm.spatial.smooth.dtype = 0;
matlabbatch{1}.spm.spatial.smooth.im = 0;
matlabbatch{1}.spm.spatial.smooth.prefix = 's';

author: translationalneuromodeling

repo: tapas


load uav3
abs = states.data;
t = 0;
n = size(abs,1);
A = [];
for i = 1:n
a1 = [t abs(i,1)  abs(i,2) abs(i,3) abs(i,4) abs(i,8) abs(i,7) abs(i,16)];
A = [A;a1];

t = t+0.01;
end

n1 = n;
plot3(A(1:n1,2), A(1:n1,3),A(1:n1,4));
axis equal


text:
function [n,msg1,msg2] = firchk(n,Fend,a,exception)
%FIRCHK   Check if specified filter order is valid.
%   FIRCHK(N,Fend,A) checks if the specified order N is valid given the
%   final frequency point Fend and the desired magnitude response vector A.
%   Type 2 linear phase FIR filters (symmetric, odd order) must have a
%   desired magnitude response vector that ends in zero if Fend = 1.  This
%   is because type 2 filters necessarily have a zero at w = pi.
%
%   If the order is not valid, a warning is given and the order
%   of the filter is incremented by one.
%
%   If A is a scalar (as when called from fircls1), A = 0 is
%   interpreted as lowpass and A = 1 is interpreted as highpass.
%
%   FIRCHK(N,Fend,A,EXCEPTION) will not warn or increase the order
%   if EXCEPTION = 1.  Examples of EXCEPTIONS are type 4 filters
%   (such as differentiators or hilbert transformers) or non-linear
%   phase filters (such as minimum and maximum phase filters).

%   Author : R. Losada
%   Copyright 1988-2004 The MathWorks, Inc.
%   $Revision: 1.7.4.5 $  $Date: 2007/12/14 15:15:06 $

matlab_v = version('-release');
matlab_v = str2double(matlab_v(1:4));

if matlab_v > 2012
    narginchk(3,4)
else
    error(nargchk(3,4,nargin,'struct'))
end

if nargin == 3,
    exception = false;
end

msg1 = '';
msg2 = '';
oddord = false; % Flag, initially we assume even order

if isempty(n) || length(n) > 1 || ~isnumeric(n) || ~isreal(n) || n~=round(n) || n<=0,
    msg1 = 'Filter order must be a real, positive integer.';
    return
end

if rem(n,2) == 1,
    oddord = true; % Overwrite flag
end

if (a(end) ~= 0) && Fend == 1 && oddord && ~exception,
    str = ['Odd order symmetric FIR filters must have a gain of zero \n'...
     'at the Nyquist frequency. The order is being increased by one.'];
    msg2 = sprintf(str);
    n = n+1;
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
% Under construction...
%
%
% Usage
% >> waverage = medianbin(ERP, bop)
%
% For Bin Operations GUI/scripting purposes
%
% Author: Javier
%

function  medb = medianbin(ERP, bop)

if nargin<1
      help medianbin
      return
end
if nargin<2
      error('ERPLAB says: You must specify bin indexes to be processed!')
end
if ~iserpstruct(ERP)
      error('ERPLAB says: Your data structure is not an ERP structure!')
end
if ischar(bop)
      bop = str2num(char(regexp(bop,'\d+','match')'))';
end

binarray = unique_bc2(bop);

if length(binarray)~=length(bop)
      fprintf('\n*** WARNING: Repeated bins were ignored.\n\n')
end
if length(binarray)<2
      error('ERPLAB says: You must specify 2 bin indexes at leat!')
end
if max(binarray)>ERP.nbin
      error('ERPLAB says: Some specified bins do not exist!')
end

datavg = ERP.bindata(:,:,binarray);
medb   = median(datavg, 3);


text:
function wathen_test05 ( )

%*****************************************************************************80
%
%% WATHEN_TEST05 measures the storage needed for the Wathen system.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 August 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WATHEN_TEST05\n' );
  fprintf ( 1, '  For various problem sizes and storage schemes, \n' );
  fprintf ( 1, '  measure the storage used for the Wathen system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                                   Predicted  Observed\n' );
  fprintf ( 1, '                              GE        Band      Band' );
  fprintf ( 1, '      Band    Sparse\n' );
  fprintf ( 1, '    NX  Elements   Nodes   storage     width     width' );
  fprintf ( 1, '   storage   storage\n' );
  fprintf ( 1, '\n' );

  nx = 1;
  ny = 1;

  for test = 1 : 6
%
%  Compute the number of unknowns.
%
    n = wathen_order ( nx, ny );
%
%  Predict the bandwidth.
%
    [ bl, bd, bu ] = wathen_bandwidth ( nx, ny );
    bw1 = bl + bd + bu;
%
%  Compute the matrix.
%
    seed = 123456789;
    [ a, seed ] = wathen_ge ( nx, ny, n, seed );
    [ na1, na2 ] = size ( a );
    storage_ge = na1 * na2;
    [ bw2, lb, db, ub ] = bandwidth ( na1, na2, a );
    storage_gb = ( 2 * lb + 1 + ub ) * n;
    nnz = length ( find ( a ~= 0.0 ) );
%
%  Report.
%
    fprintf ( 1, '  %4d      %4d  %6d  %8d  %8d  %8d  %8d  %8d\n', ...
      nx, nx * ny, n, storage_ge, bw1, bw2, storage_gb, nnz );
%
%  Ready for next iteration.
%
    nx = nx * 2;
    ny = ny * 2;

  end

  return
end


text:
function [ mW ] = CalcSpatialWeights( kernelRadius, spatialStd )
% ----------------------------------------------------------------------------------------------- %
% [ mW ] = CalcSpatialWeights( kernelRadius, spatialStd )
%   Calculates the Spatial Weights of Bilateral Filter.
% Input:
%   - vX            -   Input Vector.
%                       The point the gradient is caclaulated at.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - hObjFun       -   Objective Function.
%                       Function handler which evaluates the Objective
%                       Function at a given point - hObjFun(vX).
%                       Structure: Function Handler.
%                       Type: Function Handler.
%                       Range: NA.
%   - diffMode      -   Difference Mode.
%                       Sets the mode of operation of the Finite
%                       Differences Method - Forward, Backward, Central or
%                       Complex.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, 3, 4}.
% Output:
%   - vG            -   Gradient Vector.
%                       The numerical approximation of the gradient of the
%                       Objective Function at the input point 'vX'.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  A
% Remarks:
%   1.  B
% TODO:
%   1.  C.
% Release Notes:
%   -   1.0.000     05/10/2019  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

vX = [-kernelRadius:kernelRadius];

vW = exp(-(vX .* vX) / (2 * spatialStd * spatialStd));
mW = vW.' * vW;
mW = mW / sum(mW(:));

end

author: RoyiAvital

repo: StackExchangeCodes


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function weights_new = importance_weights(particle_system, new_rho)
% IMPORTANCE_WEIGHTS: Function updates the importance weights for each
% binary model based on the geometric bridge model and the empirical
% distribution

% Extract inputs to function
models     = particle_system.models;
model_val  = particle_system.model_val;
old_wts    = particle_system.weights;
old_rho    = particle_system.rho;

% Determine the number of models and trials
[n_models, ~] = size(models);

% Declare a vector to store weights
weights_new = zeros(n_models,1);

for i=1:n_models

    % Evaluate posterior ratios: pi(M) \propto exp(rho*f(M))
    post_new_rho = new_rho*model_val(i);
    post_old_rho = old_rho*model_val(i);

    % Find the weights corresponding to each model
    weights_new(i) = old_wts(i)*exp(post_new_rho - post_old_rho);

end

% Check for NaNs and set weight to zero
nan_weight = isnan(weights_new);
weights_new(nan_weight) = 0;

% Renormalize weights
weights_new = weights_new/sum(weights_new);

end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
classdef LayerAP < dagnn.Loss
  % @author: Hakan Bilen
  % 11 step average precision
  properties
    cls_index = 1
    resetLayer = false
    gtLabels = []
    scores   = []
    ids      = []
    aps      = []
    voc07    = true % 11 step
    classNames = {}
  end

  methods
    function outputs = forward(obj, inputs, params)
      if obj.resetLayer
        obj.gtLabels = [] ;
        obj.scores   = [] ;
        obj.ids      = [] ;
        obj.aps      = [] ;
        obj.resetLayer = false ;
      end

      if numel(inputs)==2
        obj.scores = [obj.scores gather(squeeze(inputs{1}(:,:,obj.cls_index,:)))];
        obj.gtLabels = [obj.gtLabels gather(squeeze(inputs{2}(:,:,obj.cls_index,:)))];
      elseif numel(inputs)>2
        scoresCur = gather(squeeze(inputs{1}(:,:,obj.cls_index,:)));
        gtLabelsCur = gather(squeeze(inputs{2}(:,:,obj.cls_index,:)));

        idsCur = gather(squeeze(inputs{3}));

        [lia,locb] = ismember(idsCur,obj.ids);

        if any(lia)
          obj.scores = [obj.scores scoresCur(~lia,:)];
          obj.gtLabels = [obj.gtLabels gtLabelsCur(~lia,:)];
          obj.ids = [obj.ids(:) ; idsCur(~lia,:)];

          nz = find(lia);
          for i=1:numel(nz)
            obj.scores(locb(nz(i)),:) = obj.scores(locb(nz(i)),:) + ...
              scoresCur(nz(i),:);
          end
        else
          obj.scores = [obj.scores scoresCur];
          obj.gtLabels = [obj.gtLabels gtLabelsCur];
          obj.ids = [obj.ids(:) ; idsCur]';
        end
      else
        error('wrong number of inputs');
      end

      obj.aps = obj.compute_average_precision();
      obj.average = 100 * mean(obj.aps);
      outputs{1} =  100 * mean(obj.aps);
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs = cell(1,numel(inputs));
      derInputs{1} = derOutputs{1} ;
      derParams = {} ;
    end

    function reset(obj)
      obj.resetLayer = true ;
%       obj.average = 0 ;
%       obj.aps = 0 ;
%       obj.gtLabels = [];
%       obj.scores   = [];
%       obj.ids      = [];
    end

    function printAP(obj)
      if isempty(obj.classNames)
        for i=1:numel(obj.aps)
          fprintf('class-%d %.1f\n',i,100*obj.aps(i)) ;
        end
      else
        for i=1:numel(obj.aps)
          fprintf('%-50s %.1f\n',obj.classNames{i},100*obj.aps(i)) ;
        end
      end
    end

    function aps = compute_average_precision(obj)
      assert(all(size(obj.scores)==size(obj.gtLabels)));
      % nImg = size(obj.scores,1);
      nCls = numel(obj.cls_index);

      aps = zeros(1,nCls);

      for c=1:nCls
        gt = obj.gtLabels(c,:);
        conf = obj.scores(c,:) ;
        if sum(gt>0)==0, continue ; end

        % compute average precision
        if obj.voc07
          [rec,prec,ap]=obj.VOC07ap(conf,gt) ;
        else
          [rec,prec,ap]=obj.THUMOSeventclspr(conf,gt) ;
        end
        aps(c) = ap;
      end
    end

    function [rec,prec,ap]=VOC07ap(obj,conf,gt)
      [~,si]=sort(-conf);
      tp=gt(si)>0;
      fp=gt(si)<0;

      fp=cumsum(fp);
      tp=cumsum(tp);

      rec=tp/sum(gt>0);
      prec=tp./(fp+tp);
      ap=0;
      for t=0:0.1:1
        p=max(prec(rec>=t));
        if isempty(p)
          p=0;
        end
        ap=ap+p/11;
      end
    end

    function [rec,prec,ap]=THUMOSeventclspr(obj,conf,gt)
      [so,sortind]=sort(-conf);
      tp=gt(sortind)==1;
      fp=gt(sortind)~=1;
      npos=length(find(gt==1));

      % compute precision/recall
      fp=cumsum(fp);
      tp=cumsum(tp);
      rec=tp/npos;
      prec=tp./(fp+tp);

      % compute average precision

      ap=0;
      tmp=gt(sortind)==1;
      for i=1:length(conf)
        if tmp(i)==1
          ap=ap+prec(i);
        end
      end
      ap=ap/npos;
    end

    function obj = LayerAP(varargin)
      obj.load(varargin) ;
      obj.loss = 'average_precision' ;
    end
  end
end

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function Lden = symbcholden(L,dense,DAt)
% Lden = symbcholden(L,dense,DAt)
%
% SYMBCHOLDEN  Creates Lden.{LAD, perm,dz, sign, first}
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi, dpr1fact

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 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, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

% ------------------------------------------------------------
% Symbolic forward Cholesky of dense columns, in order
% [LP, Q-blk, Q-norm, Q-tr]
% ------------------------------------------------------------
i1 = dense.l + 1;
i2 = i1 + length(dense.q);
LAD = [symbfwblk(L,dense.A(:,1:i1-1)), symbfwblk(L,DAt.denq),...
    symbfwblk(L,dense.A(:,i2:end)),symbfwblk(L,dense.A(:,i1:i2-1))];
% ------------------------------------------------------------
% Incremental ordering heuristic, excluding the Lorentz-trace cols
% ------------------------------------------------------------
[perm, dz] = incorder(LAD(:,1:length(dense.cols)));
% ------------------------------------------------------------
% Insert the trace cols with a "-1"-factor just after corresponding
% Lorentz-block columns
% ------------------------------------------------------------
Lden = finsymbden(LAD,perm,dz,i1);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [rVectECEF, vVectECEF, REci2Ecef] = getFixedFrameVectFromInertialVect(ut, rVectECI, bodyInfo, varargin)
% %getFixedFrameVectFromInertialVect Summary of this function goes here
% %   Detailed explanation goes here
%
    inputs = bodyInfo.getFixedFrameFromInertialFrameInputsCache();
    [rVectECEF, vVectECEF, REci2Ecef] = getFixedFrameVectFromInertialVect_alg(ut, rVectECI, inputs{:}, varargin{:});
end

% %getFixedFrameVectFromInertialVect Summary of this function goes here
% %   Detailed explanation goes here
%
%     if(~isempty(varargin))
%         vVectECI = varargin{1};
%     else
%         vVectECI = [NaN;NaN;NaN];
%     end
%
%     numElems = length(ut);
%
%     spinAngle = getBodySpinAngle(bodyInfo, ut);
%
%     cSA = reshape(cos(spinAngle),1,1,numElems);
%     sSA = reshape(sin(spinAngle),1,1,numElems);
%     zero = zeros(1,1,numElems);
%     one = zero + 1;
%
% %     R = [cos(spinAngle) -sin(spinAngle) 0;
% %          sin(spinAngle) cos(spinAngle) 0;
% %          0 0 1];
%
%     R = [cSA, -sSA, zero;
%          sSA, cSA, zero;
%          zero, zero, one];
%
% %     rVectECI = reshape(rVectECI,3,1);
%     rVectECI = reshape(rVectECI,3,1,numElems);
%
%     REci2Ecef = permute(R,[2,1,3]); %ND transpose
% %     rVectECEF = REci2Ecef * rVectECI;
%     rVectECEF = mtimesx(REci2Ecef, rVectECI);
%     rVectECEF = reshape(rVectECEF,3,numElems);
%
%     if(~any(isnan(vVectECI)))
%         rotRateRadSec = 2*pi/bodyInfo.rotperiod;
%         omegaRI = repmat([0;0;rotRateRadSec],1,1,numElems);
%         vVectECI = reshape(vVectECI,3,1,numElems);
%
% %         vVectECEF = REci2Ecef*(vVectECI - cross(omegaRI, rVectECI));
%         vVectECEF = mtimesx(REci2Ecef, (vVectECI - cross(omegaRI, rVectECI)));
%         vVectECEF = reshape(vVectECEF,3,numElems);
%     else
%         vVectECEF = repmat([NaN;NaN;NaN],1,numElems);
%     end


text:
function tests = test_imEuler3d(varargin)
%TESTIMEULER3D  One-line description here, please.
%
%   output = testImEuler3d(input)
%
%   Example
%   testImEuler3d
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

tests = functiontests(localfunctions);

function test_ball(testCase)

% create a simple ball
img = discreteBall(1:20, 1:20, 1:20, [10 10 10 6]);

% check EPC=1 for all adjacencies
epcTh = 1;
assertEqual(testCase, epcTh, imEuler3d(img));
assertEqual(testCase, epcTh, imEuler3d(img, 6));
assertEqual(testCase, epcTh, imEuler3d(img, 26));

% add a holes in the ball -> EPC=2
img = img & ~discreteBall(1:20, 1:20, 1:20, [10 10 10 3]);

% check EPC=1 for all adjacencies
epcTh  = 2;
assertEqual(testCase, epcTh, imEuler3d(img));
assertEqual(testCase, epcTh, imEuler3d(img, 6));
assertEqual(testCase, epcTh, imEuler3d(img, 26));

function test_Torus(testCase)

% create a torus, EPC=0
img = discreteTorus(1:60, 1:60, 1:60, [30 30 30 20 5 60 45]);

% check EPC=1 for all adjacencies
epcTh = 0;
assertEqual(testCase, epcTh, imEuler3d(img));
assertEqual(testCase, epcTh, imEuler3d(img, 6));
assertEqual(testCase, epcTh, imEuler3d(img, 26));

function test_3d_cube_C6(testCase)

% create a simple cube
img = false([5 5 5]);
img(2:4, 2:4, 2:4) = true;

% check EPC=1
epcTh = 1;
assertEqual(testCase, epcTh, imEuler3d(img, 6));

% add a hole in the ball -> EPC=2
img(3, 3, 3) = false;

% check EPC=2
epcTh  = 2;
assertEqual(testCase, epcTh, imEuler3d(img, 6));

function test_3d_cube_C26(testCase)

% create a simple cube
img = false([5 5 5]);
img(2:4, 2:4, 2:4) = true;

% check EPC=1
epcTh = 1;
assertEqual(testCase, epcTh, imEuler3d(img, 26));

% add a hole in the ball -> EPC=2
img(3, 3, 3) = false;

% check EPC=2
epcTh  = 2;
assertEqual(testCase, epcTh, imEuler3d(img, 26));

function test_3d_cubeRing_C6(testCase)

% create a simple cubic ring
img = false([5 5 5]);
img(2:4, 2:4, 2:4) = true;
img(3, 3, 2:4) = false;

% check EPC=0
epcTh = 0;
assertEqual(testCase, epcTh, imEuler3d(img, 6));

function test_3d_cubeRing_C26(testCase)

% create a cubic ring
img = false([5 5 5]);
img(2:4, 2:4, 2:4) = true;
img(3, 3, 2:4) = false;

% check EPC=0
epcTh = 0;
assertEqual(testCase, epcTh, imEuler3d(img, 26));

function test_3D_cubeDiagonals_C6(testCase)

% create a small 3D image with points along cube diagonal
img = false([5 5 5]);
for i = 1:5
    img(i, i, i) = true;
    img(6-i, i, i) = true;
    img(i, 6-i, i) = true;
    img(6-i, 6-i, i) = true;
end

epc6 = imEuler3d(img, 6);

epc6Th = 17;
assertEqual(testCase, epc6, epc6Th);

function test_3D_cubeDiagonals_C26(testCase)

% create a small 3D image with points along cube diagonal
img = false([7 7 7]);
for i = 2:6
    img(i, i, i) = true;
    img(8-i, i, i) = true;
    img(i, 8-i, i) = true;
    img(8-i, 8-i, i) = true;
end

epc26 = imEuler3d(img, 26);

epc26Th = 1;
assertEqual(testCase, epc26, epc26Th);

function test_3D_cubeDiagonals_touchingBorder_C26(testCase)

% create a small 3D image with points along cube diagonal
img = false([5 5 5]);
for i = 1:5
    img(i, i, i) = true;
    img(6-i, i, i) = true;
    img(i, 6-i, i) = true;
    img(6-i, 6-i, i) = true;
end

epc26 = imEuler3d(img, 26);

epc26Th = 1;
assertEqual(testCase, epc26, epc26Th);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function p07_story ( )

%*****************************************************************************80
%
%% P07_STORY prints the "story" for problem p07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This data is a simple symmetric set of 4 points,\n' );
  fprintf ( 1, '  for which it is interesting to develop the Shepard\n' );
  fprintf ( 1, '  interpolants for varying values of the exponent p.\n' );

  return
end


text:
function bc = specific_bc(xbd,ybd)
%ref_bc   Reference problem 3.2  boundary condition
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%   definition of Neil Madden test problem given in Roos book p.197.
%   IFISS function: DJS; 5 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
nobd=length(xbd);
bc=zeros(size(xbd));
k=find(ybd==-1); bc(k)=1;
k=find( xbd ==-1); bc(k)=(1-(1+ybd(k))/2).^3;
k=find( xbd ==1 ); bc(k)=(1-(1+ybd(k))/2).^2;
return


text:
function EEG = hlp_icaact(EEG)
% Recompute ica activations
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA

EEG.icaact = EEG.icaweights(:,EEG.icachansind)*EEG.icasphere*EEG.data(EEG.icachansind,:);
EEG.icaact = reshape(EEG.icaact, size(EEG.icaact,1), EEG.pnts, EEG.trials);

repo: REST


text:
%SiftPointFeature.support Support region of feature
%
% OUT = F.support(IM, W) is an image of the support region of the
% feature F, extracted from the image IM in which the feature appears.
% The support region is scaled to WxW and rotated so that the feature's
% orientation axis is upward.
%
% OUT = F.support(IMAGES, W) as above but if the features were extracted
% from an image sequence IMAGES then the feature is extracted from the
% appropriate image in the same sequence.
%
% [OUT,T] = F.support(IMAGES, W) as above but returns the pose of the feature
% as a 3x3 homogeneous transform in SE(2) that comprises the feature position
% and orientation.
%
% F.support(IM, W) as above but the support region is displayed.
%
% See also SiftPointFeature.

function [out,TT] = support(sf, images, N)

    if nargin < 3
        N = 50;
    end

    im = images(:,:,sf.image_id_);

    d = 2*sf.scale_;

    [Uo,Vo] = imeshgrid(N, N);

    T = se2(sf.u_, sf.v_, sf.theta_) * diag([d/N,d/N,1]) * se2(-N/2, -N/2);

    UV = transformp(T, [Uo(:) Vo(:)]');
    U = reshape(UV(1,:), size(Uo));
    V = reshape(UV(2,:), size(Vo));

    [Ui,Vi] = imeshgrid(im);

    im2 = interp2(Ui, Vi, idouble(im), U, V);

    if nargout == 0
        idisp(im2)
    elseif nargout == 1
        out = im2;
    elseif nargout == 2
        out = im2;
        TT = T;
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function infsup(c)
%INFSUP       Display of interval hessians in infsup notation
%
%   infsup(c)
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 02/11/06     S.M. Rump  SparseInfNanFlag removed
% modified 06/04/09     S.M. Rump  Comment
% modified 08/26/12     S.M. Rump  global variables removed
%

  if nargin<2
    name = inputname(1);
    if isempty(name)                    % happens for display(hessianinit(random))
      name = 'ans';
    end
  end

  if isa(c.x,'intval')
    display_gen(c,name,'infsup')
  else
    display_gen(c,name,'display_')
  end


text:
function []  = imwriteBW(x,filename)
	m = min(x(:));
	M = max(x(:));

	xn =  (x-m)/(M-m);
	imwrite(xn,filename);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [superpixels, clusters, Z] = doAmodalCompletion(imName, paths, ucm2, pc, amodalParam)
	ucmThresh = amodalParam.ucmThresh;
	thresh = amodalParam.thresh;

	superpixels = bwlabel(ucm2 < amodalParam.ucmThresh);
	superpixels = superpixels(2:2:end, 2:2:end);
	[clusters, Z] = amodalCompletion(superpixels, pc, thresh);

	fileName = fullfile(paths.amodalDir, strcat(imName, '.mat'));
	save(fileName, 'ucmThresh', 'superpixels', 'thresh', 'clusters', 'Z');
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function CPT = CPD_to_CPT(CPD)
% CPD_TO_CPT Convert the tabular_decision_node to a CPT
% CPT = CPD_to_CPT(CPD)

CPT = CPD.CPT;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function rad=dms2rad(dms)
% DMS2RAD  Converts degrees-minutes-seconds to radians.
%   Vectorized.
% Version: 12 Mar 00
% Useage:  rad=dms2rad(dms)
% Input:   dms - [d m s] array of angles in deg-min-sec, where
%                d = vector of degrees
%                m = vector of minutes
%                s = vector of seconds
% Output: rad - vector of angles in radians

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

d=dms(:,1);
m=dms(:,2);
s=dms(:,3);
dec=abs(d)+abs(m)./60+abs(s)./3600;
rad=dec.*pi./180;
ind=(d<0 | m<0 | s<0);
rad(ind)=-rad(ind);


text:
function [max_val,ind] = mDmax(x)
%MDMAX Largest component in a multidimensional matrix.
%   For a matrix X, MDMAX(X) is the largest element in X.
%
%   [Y,I] = MDMAX(X) returns the indices of the maximum value in vector I.
%
%   When X is complex, the maximum is computed using the magnitude
%   MAX(ABS(X)). In the case of equal magnitude elements, then the phase
%   angle MAX(ANGLE(X)) is used.
%
%   NaN's are ignored when computing the maximum. When all elements in X
%   are NaN's, then the first one is returned as the maximum.
%
%   Example: If X = [2 8 4;   then mDmax(X) is 9,
%                    7 3 9]
%   and
%            [val,ind] = mDmax(X) returns val=9 and ind=[2 3]
%
%   Note: The vector ind can be used to access elements from X using this
%   neat trick:
%              sub = num2cell(ind)
%              val = x(sub{:});

[max_val,position] = max(x(:));
ind = myind2sub(size(x),position);

%--------------------------------------------------------------
function ind = myind2sub(siz,ndx)

dim = numel(siz);
sub = cell(dim,1);
[sub{:}] = ind2sub(siz,ndx);
ind = cell2mat(sub);

% Old implementation
% siz = double(siz);
%
% n = length(siz);
% k = [1 cumprod(siz(1:end-1))];
% for i = n:-1:1,
%     vi = rem(ndx-1, k(i)) + 1;
%     vj = (ndx - vi)/k(i) + 1;
%     ind(:,i) = vj;
%     ndx = vi;
% end


text:
files = {
    fullfile('VPibv_10_11_02','calibration_CenterSpellerMVEP_VPibv')
    fullfile('VPibq_10_09_24','calibration_CenterSpellerMVEP_VPibq')
    fullfile('VPiac_10_10_13','calibration_CenterSpellerMVEP_VPiac')
    };

%% Define some settings
disp_ival= [-200 1000];
ref_ival= [-200 0];
crit_maxmin= 70;
crit_ival= [100 800];
crit_clab= {'F9,z,10','AF3,4'};
clab= {'Cz','PO7'};
colOrder= [1 0 1; 0.4 0.4 0.4];

%% load all data sets
nsub = length(files);
for isub = 1:nsub
    file = files{isub};
    file= fullfile(BTB.DataDir, 'demoMat', file);

    %% Load data
    fprintf('loading %s \n', file)

    [cnt, mrk, mnt] = file_loadMatlab(file);

    %% Apply highpass filter to reduce drifts
    b= procutil_firlsFilter(0.5, cnt.fs);
    cnt= proc_filtfilt(cnt, b);

    %% Artifact rejection based on variance criterion
    mrk= reject_varEventsAndChannels(cnt, mrk, disp_ival, 'verbose', 1);

    %% Segmentation
    epo= proc_segmentation(cnt, mrk, disp_ival);

    %% Artifact rejection based on maxmin difference criterion on frontal chans
    [epo iArte] = proc_rejectArtifactsMaxMin(epo, crit_maxmin, 'Clab',crit_clab, ...
        'Ival',crit_ival, 'Verbose',1);

    %% Baseline subtraction, and calculation of a measure of discriminability
    epo= proc_baseline(epo, ref_ival);

    epos_av{isub} = proc_average(epo, 'Stats', 1);

    % three different but almost equivalent ways to make statistics about class differences
    epos_diff{isub} = proc_classmeanDiff(epo, 'Stats', 1);
    epos_r{isub} = proc_rSquareSigned(epo, 'Stats', 1);
    epos_auc{isub} = proc_aucValues(epo, 'Stats', 1);

end

%% grand average
epo_av = proc_grandAverage(epos_av, 'Average', 'INVVARweighted', 'Stats', 1, 'Bonferroni', 1, 'Alphalevel', 0.01);
epo_r = proc_grandAverage(epos_r, 'Average', 'INVVARweighted', 'Stats', 1, 'Bonferroni', 1, 'Alphalevel', 0.01);
epo_diff = proc_grandAverage(epos_diff, 'Average', 'INVVARweighted', 'Stats', 1, 'Bonferroni', 1, 'Alphalevel', 0.01);
epo_auc = proc_grandAverage(epos_auc, 'Average', 'INVVARweighted', 'Stats', 1, 'Bonferroni', 1, 'Alphalevel', 0.01);

%% Select some discriminative intervals, with constraints to find N2, P2, P3 like components.
constraint= ...
    {{-1, [100 300], {'I#','O#','PO7,8','P9,10'}, [50 300]}, ...
    {1, [200 350], {'P3-4','CP3-4','C3-4'}, [200 400]}, ...
    {1, [400 500], {'P3-4','CP3-4','C3-4'}, [350 600]}};
[ival_scalps, nfo]= ...
    procutil_selectTimeIntervals(epo_r, 'Visualize', 0, 'VisuScalps', 1, ...
    'Title', util_untex(file), ...
    'Clab',{'not','E*'}, ...
    'Constraint', constraint);
%printFigure('r_matrix', [18 13]);
ival_scalps= visutil_correctIvalsForDisplay(ival_scalps, 'fs',epo.fs);

%% plot classwise grand-average ERPs
fig_set(1);
H= plot_scalpEvolutionPlusChannel(epo_av, mnt, clab, ival_scalps, defopt_scalp_erp, ...
    'ColorOrder',colOrder);
grid_addBars(epo_r);
%printFigure(['erp_topo'], [20  4+5*size(epo.y,1)]);

%% plot difference of the class means
fig_set(2, 'Resize',[1 2/3]);
plot_scalpEvolutionPlusChannel(epo_diff, mnt, clab, ival_scalps, defopt_scalp_r);
%printFigure(['erp_topo_r'], [20 9]);

% plot signed log10 p-values of the null hypothesis
% that the difference of the class means is zero
% interpretation: abs(sgnlogp) > 1   <-->  p < 0.1
%                 abs(sgnlogp) > 2   <-->  p < 0.01
%                 abs(sgnlogp) > 3   <-->  p < 0.001 , and so on
fig_set(3, 'Resize',[1 2/3]);
epo_diff_sgnlogp = epo_diff;
epo_diff_sgnlogp.x = epo_diff_sgnlogp.sgnlogp;
epo_diff_sgnlogp.yUnit = 'sgnlogp';
plot_scalpEvolutionPlusChannel(epo_diff_sgnlogp, mnt, clab, ival_scalps, defopt_scalp_r);
%printFigure(['erp_topo_r'], [20 9]);

%% now plot differences again, with all insignificant results set to zero
epo_diff.x = epo_diff.x.*epo_diff.sigmask;
fig_set(4, 'Resize',[1 2/3]);
plot_scalpEvolutionPlusChannel(epo_diff, mnt, clab, ival_scalps, defopt_scalp_r);

author: bbci


text:
function amps = mv_amps(mv,runs);
%
% amps = mv_amps(mv,[runs]);
%
% Return a 2D matrix of format
% voxels x conditions containing
% response amplitudes, evaluated
% according to the 'ampType' parameter
% set in the params struct. Does not compute amps
% for the null (0) condition.
%
% ras, 05/05.

if notDefined('mv')
    mv = get(gcf,'UserData');
end

allRuns = unique(mv.trials.run);

if notDefined('runs')
    runs = allRuns;
end

% check that we have an amplitude type defined
if ~checkfields(mv,'params','ampType')
    % set up a dialog, get it
    ui.string = 'Method to Calculate Amplitudes?';
    ui.fieldName = 'ampType';
    ui.list = {'Peak-Bsl Difference', 'GLM Betas', ...
        'Dot-product Relative Amps' 'raw'};
    ui.style = 'popup';
    ui.value = ui.list{1};

    resp = generalDialog(ui,'Select Amplitude Type');
    ampInd = cellfind(ui.list,resp.ampType);
    opts = {'difference' 'betas' 'relamps' 'raw'};
    mv.params.ampType = opts{ampInd};
end

nConds = sum(mv.trials.condNums > 0);

switch mv.params.ampType
    case 'difference',
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % peak-bsl difference              %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        [tSeries, trials] = er_dataSubset(mv, runs);
        voxData = er_voxDataMatrix(tSeries, trials, mv.params);

        voxAmps = er_voxAmpsMatrix(voxData, mv.params);
        amps = permute(nanmeanDims(voxAmps,1), [2 3 1]);

    case 'betas',
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % apply a GLM and get beta values  %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Get data from selected runs only
        [tSeries, trials] = er_dataSubset(mv, runs);

        % Build a design matrix, apply the glm, grab betas
        [X, nh, hrf] = glm_createDesMtx(trials, mv.params, tSeries, 0);

        model = glm(double(tSeries), X, mv.params.framePeriod, nh);
        amps = permute(model.betas(1,1:nConds,:), [3 2 1]);

%         % also add dc components to each beta value, to ensure the shift is
%         % not lost
%         dc = permute( model.betas(1,nConds+1:end,:), [3 2 1]);
%         amps = amps + repmat(mean(dc,2), [1 nConds]);

    case 'trialbetas'
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % apply a GLM and get beta values for each trial %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [tSeries, trials] = er_dataSubset(mv, runs);
        nConds = sum(trials.cond > 0);
        counter = 0;
        for i = 1:length(trials.cond)
            if (trials.cond(i))
                counter = counter + 1;
                trials.cond(i) = counter;
            end
        end

        [X, nh, hrf] = glm_createDesMtx(trials, mv.params, tSeries, 0);

        model = glm(double(tSeries), X, mv.params.framePeriod, nh);
        amps = permute(model.betas(1,1:nConds,:), [3 2 1]);
        return;

    case 'raw',
		[tSeries, trials] = er_dataSubset(mv, runs);
		mv.params.normBsl=0;
		mv.params.ampType='raw'
		mv.voxData = er_voxDataMatrix(tSeries,mv.trials,mv.params);
		mv.mvAmps=er_voxAmpsMatrix(mv.voxData, mv.params);
		amps = permute(nanmeanDims(mv.voxAmps(runs,:,:),1), [2 3 1]);
    case 'zscore',
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Mean peak-bsl / std. deviation   %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        tmp = mv;
        tmp.params.ampType = 'difference';
        mu = mv_amps(tmp, runs);
        sigma = mv_stdev(tmp, runs);
        amps = mu ./ sigma;
        return      % already chose selected conditions

    case 'relamps',
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % dot-product relative amplitudes: %
        % * not currently implemented *    %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        amps = er_relamps(mv.voxData(:,runs,:,:));
        amps = squeeze(nanmean(amps));

        % return amplitudes for selected conditions only
        sel = find(tc_selectedConds(mv));
        amps = amps(:,sel-1);

    case 'deconvolved',
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Amplitudes of deconvolved time courses %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Get data from selected runs only
        [tSeries, trials] = er_dataSubset(mv, runs);

        % set up and run a GLM
        params = mv.params;
        params.glmHRF = 0;
        TR = mv.params.framePeriod;
        [X nh] = glm_createDesMtx(trials, params, tSeries);
        model = glm(tSeries, X, TR, nh);

        % get peak-bsl amplitudes from the deconvolved time courses
        frameWindow = unique(round(mv.params.timeWindow./TR));
        prestim = -1 * frameWindow(1);
        peakFrames = unique(round(mv.params.peakPeriod./TR));
        bslFrames = unique(round(mv.params.bslPeriod./TR));
        peakFrames = find(ismember(frameWindow,peakFrames));
        bslFrames = find(ismember(frameWindow,bslFrames));

        nConds = size(model.betas, 2);

        if mv.params.normBsl==1
            offset = mean(model.betas(bslFrames,:,:), 1);
            offset = repmat(offset, [length(frameWindow) 1 1]);
            model.betas = model.betas - offset;
        end

        for i = 1:nConds
            bsl = model.betas(bslFrames, i, :);
            peak = model.betas(peakFrames, i, :);
            amps(:,i) = squeeze(nanmean(peak) - nanmean(bsl))';
		end
	otherwise
		fprintf(1,'error ampType %s does not exist \n', mv.params.ampType);
end

% return amplitudes for selected conditions only
sel = find(tc_selectedConds(mv));
amps = amps(:,sel-1);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcMesh(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

figure;
mesh(double(handles.BMatrix));
colormap(handles.V.Color_map);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [FluxCorrelations, PValues, TaxonomyInfo] = correlateFluxWithTaxonAbundance(abundancePath, fluxPath, infoFilePath, corrMethod)
% Part of the Microbiome Modeling Toolbox. This function calculates and
% plots the correlations between fluxes for one or more reactions of
% interest in a number of microbiome samples and the relative microbe
% abundance on different taxonomical levels in the same samples.
% The function should be used after running mgPipe to identify correlations
% between the computed metabolic profiles and specific taxa in the samples.
%
% USAGE
%
%     [FluxCorrelations, PValues, TaxonomyInfo] = correlateFluxWithTaxonAbundance(abundancePath, fluxPath, taxonomy, corrMethod)
%
% INPUTS:
%    abundancePath:     Path to the .csv file with the abundance data.
%                       Needs to be in same format as example file
%                       'cobratoolbox/papers/018_microbiomeModelingToolbox/examples/normCoverage.csv'
%     fluxPath:         Path to the .csv file with the fluxes for reactions
%                       of interest with sample IDs as rows and reaction
%                       IDs in microbiome community models as columns
%
% OPTIONAL INPUTS:
%    infoFilePath:      Path to the spreadsheet with the taxonomy information
%                       on organisms (default: AGORA_infoFile.xlsx)
%     corrMethod:       Method to compute the linear correlation
%                       coefficient. Allowed inputs: 'Pearson' (default),
%                       'Kendall', 'Spearman'.
%
% OUTPUTS:
%     FluxCorrelations: Structure with correlations between fluxes for each
%                       reaction and abundances on taxon levels
%     PValues:          p-values corresponding to each calculated
%                       correlation
%     TaxonomyInfo:     Taxonomical information on each taxon level
%
% .. Author: Almut Heinken, 03/2018
%                           10/2018:  changed input to location of the csv file with the
%                                     abundance data
%                           01/2020:  adapted to be suitable for pan-models, and
%                                     changed flux input to a csv file.

% read the csv file with the abundance data
abundance = readInputTableForPipeline(abundancePath);
if isnumeric(abundance{2, 1})
    abundance(:, 1) = [];
end
% adapt IDs if neccessary
abundance(1,2:end) = strrep(abundance(1,2:end),'-','_');

fluxes = readInputTableForPipeline(fluxPath);
fluxes(1,:) = strrep(fluxes(1,:),'microbiota_model_diet_','');
fluxes(1,:) = strrep(fluxes(1,:),'microbiota_model_samp_','');
fluxes(:,1) = strrep(fluxes(:,1),'microbiota_model_diet_','');
fluxes(:,1) = strrep(fluxes(:,1),'microbiota_model_samp_','');

% check if data is from same samples
if ~isempty(setdiff(fluxes(1,2:end),abundance(1,2:end)))
    fluxes=fluxes';
    % if it still does not match
    if ~isempty(setdiff(fluxes(1,2:end),abundance(1,2:end)))
        error('Sample IDs in abundance and flux files do not agree!')
    end
end

% load database
database=loadVMHDatabase;

fluxes(:,1)=strrep(fluxes(:,1),'EX_','');
fluxes(:,1)=strrep(fluxes(:,1),'(e)','');
fluxes(:,1)=strrep(fluxes(:,1),'[fe]','');
% for i=2:size(fluxes,1)
% fluxes{i,1}=database.metabolites{find(strcmp(database.metabolites(:,1),fluxes{i,1})),2};
% end

% Get the taxonomy information
if exist('infoFilePath','var')
    taxonomy = readInputTableForPipeline(infoFilePath);
else
    taxonomy = readInputTableForPipeline('AGORA_infoFile.xlsx');
end

if ~exist('corrMethod', 'var')  % Define correlation coefficient method if not entered
    corrMethod = 'Pearson';
end

% Calculate the abundance in each sample on all taxon levels
TaxonomyLevels = {
    'Phylum'
    'Class'
    'Order'
    'Family'
    'Genus'
    'Species'
};
% extract the list of entries on each taxonomical level and prepare the
% summarized abundance table
fprintf('Calculating the relative abundances on all taxon levels. \n')
for t = 1:size(TaxonomyLevels, 1)
    % find the columns corresponding to each taxonomy level and the list of
    % unique taxa
    taxonCol = find(strcmp(taxonomy(1, :), TaxonomyLevels{t}));
    % find and save all entries
    taxa = unique(taxonomy(2:end, taxonCol));
    % exclude unclassified entries
    taxa(strncmp('unclassified', taxa, taxonCol)) = [];
    TaxonomyLevels{t, 2} = taxa;
    for i = 1:length(taxa)
        SampleAbundance.(TaxonomyLevels{t}){1, i + 1} = taxa{i};
        for j = 2:size(abundance, 2)
            SampleAbundance.(TaxonomyLevels{t}){j, 1} = abundance{1, j};
            SampleAbundance.(TaxonomyLevels{t}){j, i + 1} = 0;
        end
    end
end

% Go through the abundance data and summarize taxon abundances for each
% strain in at least one sample

% Find the right column for the input data (strains, species,..)
if sum(strncmp(abundance(:,1),'pan',3))==size(abundance,1)-1
abundance(:,1)=regexprep(abundance(:,1),'pan','','once');
end
inputTaxa={};
for i=1:size(taxonomy,2)
    if ~isnumeric(taxonomy{2,i})
        taxa=strrep(taxonomy(:,i),' ','_');
        taxa=strrep(taxa,'.','_');
        taxa=strrep(taxa,'/','_');
        taxa=strrep(taxa,'-','_');
        taxa=strrep(taxa,'-','_');
        taxa=strrep(taxa,'[','');
        taxa=strrep(taxa,']','');
        taxa=strrep(taxa,')','');
        taxa=strrep(taxa,'__','_');
        if length(intersect(abundance(2:end,1),taxa))==size(abundance,1)-1
            inputTaxa=taxa;
            inputCol=i;
            break
        end
    end
end
if isempty(inputTaxa)
    error('Some taxa in the abundance file are not found in the taxonomy file!')
end

for i = 2:size(abundance, 2)
    for j = 2:size(abundance, 1)
        for t = 1:size(TaxonomyLevels, 1)
            % find the taxon for the current strain
            taxonCol = find(strcmp(taxonomy(1, :), TaxonomyLevels{t}));
            if taxonCol >= inputCol
                findTax = taxonomy(find(strcmp(abundance{j, 1}, inputTaxa)), taxonCol);
                if isempty(strfind(findTax{1}, 'unclassified'))
                    % find the taxon for the current strain in the sample abundance
                    % variable
                    findinSampleAbun = find(strcmp(findTax{1}, SampleAbundance.(TaxonomyLevels{t})(1, :)));
                    % sum up the relative abundance
                    if contains(version,'(R202') % for Matlab R2020a and newer
                        SampleAbundance.(TaxonomyLevels{t}){i, findinSampleAbun} = SampleAbundance.(TaxonomyLevels{t}){i, findinSampleAbun} + abundance{j, i};
                    else
                        SampleAbundance.(TaxonomyLevels{t}){i, findinSampleAbun} = SampleAbundance.(TaxonomyLevels{t}){i, findinSampleAbun} + str2double(abundance{j, i});
                    end
                end
            end
        end
    end
end
% remove the taxa not present in samples or only present in small abundances
for t = 1:size(TaxonomyLevels, 1)
    delArray = [];
    cnt = 1;
    for i = 2:size(SampleAbundance.(TaxonomyLevels{t}), 2)
        for j = 2:size(SampleAbundance.(TaxonomyLevels{t}), 1)
            abun(j - 1, 1) = SampleAbundance.(TaxonomyLevels{t}){j, i};
        end
        if sum(abun) < 0.005
            delArray(cnt, 1) = i;
            cnt = cnt + 1;
        end
    end
    SampleAbundance.(TaxonomyLevels{t})(:, delArray) = [];
end
% find the flux data for each reaction
fprintf('Calculating the correlations between fluxes and abundances. \n')
for i = 2:size(fluxes, 1)
    data = [];
    for m = 2:size(fluxes, 2)
        data(m - 1, 1) = str2double(string(fluxes{i, m}));
    end
    for t = 1:size(TaxonomyLevels, 1)
        FluxCorrelations.(TaxonomyLevels{t}){1, i} = fluxes{i, 1};
        PValues.(TaxonomyLevels{t}){1, i} = fluxes{i, 1};
        % find the abundance data for each taxon
        for j = 2:size(SampleAbundance.(TaxonomyLevels{t}), 2)
            FluxCorrelations.(TaxonomyLevels{t}){j, 1} = SampleAbundance.(TaxonomyLevels{t}){1, j};
            PValues.(TaxonomyLevels{t}){j, 1} = SampleAbundance.(TaxonomyLevels{t}){1, j};
            % find the abundance data for each sample
            dataTaxa = data;
            for k = 2:size(SampleAbundance.(TaxonomyLevels{t}), 1)
                % match with correct individual in flux table
                sampleInFluxes = find(strcmp(fluxes(1, :), SampleAbundance.(TaxonomyLevels{t}){k, 1}));
                dataTaxa(sampleInFluxes - 1, 2) = SampleAbundance.(TaxonomyLevels{t}){k, j};
            end
            % exclude NaNs
            dataTaxa(find(isnan(dataTaxa(:,1))),:)=[];
            dataTaxa(find(isnan(dataTaxa(:,2))),:)=[];
            % calculate the correlation with the given correlation coefficient method
            [RHO, PVAL] = corr(dataTaxa(:, 1), dataTaxa(:, 2), 'type', corrMethod);
            if isnan(RHO)
                RHO = 0;
            end
            if abs(RHO) < 0.0000000001
                RHO = 0;
            end
            FluxCorrelations.(TaxonomyLevels{t}){j, i} = RHO;
            PValues.(TaxonomyLevels{t}){j, i} = PVAL;
        end
    end
end

% remove entries that are only weak correlations and empty taxonomy entries
for t = 1:size(TaxonomyLevels, 1)
    cnt=1;
    delArray=[];
    for j=2:size(FluxCorrelations.(TaxonomyLevels{t}),2)
        if ~any(abs(cell2mat(FluxCorrelations.(TaxonomyLevels{t})(2:end,j))) > 0.3)
            delArray(cnt,1)=j;
            cnt=cnt+1;
        end
    end
    FluxCorrelations.(TaxonomyLevels{t})(:,delArray)=[];

    cnt=1;
    delArray=[];
    for j=2:size(FluxCorrelations.(TaxonomyLevels{t}),1)
        if ~any(abs(cell2mat(FluxCorrelations.(TaxonomyLevels{t})(j,2:end))) > 0.3)
            delArray(cnt,1)=j;
            cnt=cnt+1;
        end
    end
    FluxCorrelations.(TaxonomyLevels{t})(delArray,:)=[];

    FluxCorrelations.(TaxonomyLevels{t})(find(strcmp(FluxCorrelations.(TaxonomyLevels{t})(:,1),'')),:)=[];
end

% translate to metabolite descriptions
for t = 2:size(TaxonomyLevels, 1)
    for i=2:size(FluxCorrelations.(TaxonomyLevels{t}),2)
        try
            FluxCorrelations.(TaxonomyLevels{t}){1,i}=database.metabolites{find(strcmp(database.metabolites(:,1),FluxCorrelations.(TaxonomyLevels{t}){1,i})),2};
        catch
            warning('Flux label could not be translated to metabolite name.')
        end
    end
end

% export taxonomical information
taxonCol = 'Phylum';
% remove unnecessary columns
taxonomy(:,taxonCol+1:end)=[];

for t = 2:size(TaxonomyLevels, 1)
    TaxonomyReduced=taxonomy;
    taxonCol = find(strcmp(taxonomy(1, :), TaxonomyLevels{t}));
    TaxonomyReduced(:,1:taxonCol-1)=[];
    % remove duplicate entries
    [C,IA] = unique(TaxonomyReduced(:,1),'stable');
    % remove unclassified taxa
    findUncl=find(contains(C,'unclassified'));
    IA(findUncl,:)=[];
    TaxonomyInfo.(TaxonomyLevels{t})=TaxonomyReduced(IA,:);
end

set(0, 'DefaultTextInterpreter', 'none')

% Plot the calculated correlations.
for t = 1:length(TaxonomyLevels)
    data = string(FluxCorrelations.(TaxonomyLevels{t})(2:end, 2:end));
    data = str2double(data);

    if size(FluxCorrelations.(TaxonomyLevels{t}),2) > 5 && size(FluxCorrelations.(TaxonomyLevels{t}),1) > 5
        if size(FluxCorrelations.(TaxonomyLevels{t}),2) < 50
            cgo = clustergram(data,...
                'RowLabels', FluxCorrelations.(TaxonomyLevels{t})(2:end,1),...
                'ColumnLabels', FluxCorrelations.(TaxonomyLevels{t})(1,2:end),...
                'ColumnLabelsRotate', 45, ...
                'Cluster', 'all', ...
                'ShowDendrogram','True', ...
                'symmetric','False', ...
                'colormap', 'turbo' ...
                );
        else
            cgo = clustergram(data,...
                'RowLabels', FluxCorrelations.(TaxonomyLevels{t})(2:end,1),...
                'Cluster', 'all', ...
                'ShowDendrogram','True', ...
                'symmetric','False', ...
                'colormap', 'turbo' ...
                );
        end
        h = plot(cgo);
        set(h,'TickLabelInterpreter','none');
        colorbar(h)
    end
end

end

author: opencobra

repo: cobratoolbox


text:
function data = metoffice_get_data(dataset, anomalies, remove_val)

if nargin < 3 || isempty(remove_val)
    remove_val = 1;
end

%
% Load the data and remove land areas
%

% $$$ logname = getenv('LOGNAME');
% $$$ if strcmp( logname, 'alexilin' ) || strcmp( logname, 'hadia' )
% $$$     jaakkopath
% $$$ end
% $$$
% $$$ os = system_dependent('getos');
% $$$ infinland = ~isempty(strfind(os,'Ubuntu')) || ~isempty(strfind(os,'Linux'));
% $$$ if infinland
% $$$     datapath = '/share/climate/data/UK_Met_Office/RecTest/DATASETS';
% $$$     % Folder for saving the results
% $$$     if strcmp( logname, 'alexilin' )
% $$$         folder = '/share/climate/data/UK_Met_Office/RecTest/RESULTS';
% $$$     else
% $$$         folder = '/share/work/jluttine/metoffice/rectest';
% $$$     end
% $$$ else
% $$$     datapath = '/net/home/h03/hadia/data/RecTest/';
% $$$     folder = '/net/home/h03/hadia/data/RecTest/RESULTS';
% $$$ end

datapath = metoffice_get_path();

switch dataset
case {'hadsst2d1','hadsst2d2','mohsst5d1','mohsst5d2'}
    datafile = [ datapath, dataset ];
    maskfile = 'mask';
case {'aari','aari8'}
    datafile = [ datapath, dataset ];
    maskfile = 'mask_seaice';
otherwise
    error(sprintf('Unknown dataset %s',dataset));
end

data = load(datafile);

[LON,LAT] = meshgrid(data.lon,data.lat);
data.coordinates = double(metoffice_remove_mask([LON(:),LAT(:)], dataset))';
%data.coordinates = double(metoffice_remove_bins([LON(:),LAT(:)],maskfile))';

%if isempty( findstr( maskfile, 'seaice' ) )

    % SST
    Itrain = ~isnan(data.sst);
    if remove_val
        Ival = metoffice_get_valset( remove_val, data.time, Itrain );
        %Ival = get_valset( remove_val, data.time, Itrain );
        data.sst(Ival) = NaN;
        fprintf( 'Validation set %d removed\n', remove_val )
    end

    % Remove land areas
    %data.data = double(metoffice_remove_bins(data.sst,maskfile));
    data.data = double(metoffice_remove_mask(data.sst, dataset));
    data = rmfield( data, 'sst' );

% $$$ else
% $$$
% $$$     % Sea ice
% $$$     % Remove areas
% $$$     data.data = double(metoffice_remove_bins(data.sice,maskfile));
% $$$     data = rmfield( data, 'sice');
% $$$
% $$$ end

data.clim = double(metoffice_remove_mask(data.clim, dataset));
data.gridsize = metoffice_remove_mask(cosd(LAT(:)), dataset);
% $$$ data.clim = double(metoffice_remove_bins(data.clim,maskfile));
% $$$ data.gridsize = metoffice_remove_bins(cosd(LAT(:)),maskfile);

if anomalies==true
  % Remove climatological monthly averages
  dv = datevec( data.time ); month = dv(:,2);

  data.data = data.data - data.clim(:,month);
  clear month dv

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% std_checkfiles() - Check all STUDY files consistency
%
% Usage:
%                >> boolval = std_checkfiles(STUDY, ALLEEG);
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - All EEGLAB datasets
%
% Outputs:
%   boolval    - [0|1] 1 if uniform
%
% Authors: Arnaud Delorme, CERCO, 2010-

% Copyright (C) Arnaud Delorme, CERCO
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [boolval npersubj] = std_checkfiles(STUDY, ALLEEG);

if nargin < 2
    help std_checkfiles;
    return;
end;
return;

filetypes = { 'daterp' 'datspec' 'datersp' 'datitc' 'dattimef' ...
              'icaerp' 'icaspec' 'icaersp' 'icaitc' 'icatopo' };

% set channel interpolation mode
% ------------------------------
uniformchannels = std_uniformsetinds( STUDY );

disp('---------------------------------------------');
disp('Checking data files integrity and consistency');
passall = 1;
for index = 1:length(filetypes)

    % scan datasets
    % -------------
    [ tmpstruct compinds filepresent ] = std_fileinfo(ALLEEG, filetypes{index});
    if ~isempty(tmpstruct)
        fprintf('Files of type "%s" detected, checking...',  filetypes{index});
    end;

    % check if the structures are equal
    % ---------------------------------
    notequal = any(~filepresent);
    if ~isempty(tmpstruct)
        if any(filepresent == 0)
            fprintf(' Error, some files inconsistent or missing\n');
            notequal = 1;
            passall = 0;
        else
            firstpass = 1;
            fields = fieldnames(tmpstruct);
            for f_ind = 1:length(fields)

                firstval = getfield(tmpstruct, {1}, fields{f_ind});
                for dat = 2:length(tmpstruct)
                    tmpval = getfield(tmpstruct, {dat}, fields{f_ind});
                    if ~isequal(firstval, tmpval)

                        % check for NaNs
                        if iscell(firstval) && iscell(tmpval)
                            for cind = 1:length(firstval)
                                if isreal(firstval{cind}) && ~isempty(firstval{cind}) && isnan(firstval{cind}(1))
                                    firstval{cind} = 'NaN';
                                end;
                            end;
                            for cind = 1:length(tmpval)
                                if isreal(tmpval{cind}) && ~isempty(tmpval{cind}) && isnan(tmpval{cind}(1))
                                    tmpval{cind} = 'NaN';
                                end;
                            end;
                        end;

                        if ~isequal(firstval, tmpval)
                            if ~strcmpi(fields{f_ind}, 'labels') || strcmpi(uniformchannels, 'on')
                                if firstpass == 1, fprintf('\n'); firstpass = 0; end;
                                fprintf('  Error, difference accross data files for field "%s"\n', fields{f_ind});
                                notequal = 1;
                                passall = 0;
                                break;
                            end;
                        end;
                    end;
                end;
            end;
        end;
    end;

    % check the consistency of changrp and cluster with saved information
    % -------------------------------------------------------------------
    if isempty(tmpstruct), notequal = 1; end;
    if filetypes{index}(1) == 'd' && notequal == 0
        % scan all channel labels
        % -----------------------
        if isfield(tmpstruct(1), 'labels')
            for cind = 1:length(STUDY.changrp)
                if notequal == 0
                    for inddat = 1:length(ALLEEG)
                        tmpind = cellfun(@(x)(find(x == inddat)), STUDY.changrp(cind).setinds(:), 'uniformoutput', false);
                        indnonempty = find(~cellfun(@isempty, tmpind(:)));
                        if ~isempty(indnonempty)
                            tmpchan = STUDY.changrp(cind).allinds{indnonempty}(tmpind{indnonempty}); % channel index for dataset inddat

                            tmpchan2 = strmatch(STUDY.changrp(cind).name, tmpstruct(inddat).labels, 'exact'); % channel index in file
                            if ~isempty(tmpchan2) || ~strcmpi(filetypes{index}, 'datspec') % the last statement is because channel labels did not use to be saved in spec files
                                if ~isequal(tmpchan2, tmpchan)
                                    fprintf('\nError: channel index in STUDY.changrp(%d) for dataset %d is "%d" but "%d" in data files\n', cind, inddat, tmpchan, tmpchan2);
                                    notequal = 1;
                                    break;
                                end;
                            end;
                        end;
                    end;
                end;
            end;
        end;
    elseif notequal == 0 && ~isempty(STUDY.cluster) % components
        % check that the cell structures are present
        % ------------------------------------------
        if ~isfield(STUDY.cluster, 'setinds')
            STUDY.cluster(1).setinds = [];
            STUDY.cluster(1).allinds = [];
        end;
        for cind = 1:length(STUDY.cluster)
            if isempty(STUDY.cluster(cind).setinds)
                STUDY.cluster(cind) = std_setcomps2cell(STUDY, cind);
            end;
        end;
        for cind = 1:length(STUDY.cluster)
            if notequal == 0
                for inddat = 1:length(ALLEEG)
                    tmpind = cellfun(@(x)(find(x == inddat)), STUDY.cluster(cind).setinds(:), 'uniformoutput', false);
                    indnonempty = find(~cellfun(@isempty, tmpind(:)));
                    tmpcomp = [];
                    for jind = 1:length(indnonempty)
                        tmpcomp = [ tmpcomp STUDY.cluster(cind).allinds{indnonempty(jind)}(tmpind{indnonempty(jind)}) ];
                    end;

                    if ~isempty(setdiff(tmpcomp, compinds{inddat}))
                        if ~(isempty(compinds{inddat}) && strcmpi(filetypes{index}, 'icatopo'))
                            fprintf('\nError: some components in clusters %d are absent from .%s files\n', cind, filetypes{index});
                            notequal = 1;
                            passall  = 0;
                            break;
                        end;
                    end;
                end;
            end;
        end;
    end;
    if notequal == 0, fprintf(' Pass\n'); end;
end;

if ~passall
    disp('**** Recompute any measure above not receiving a "Pass" by')
    disp('**** calling menu items "STUDY > Precompute Channel/Component measures" ');
    disp('**** and by selecting the "Recompute even if present on disk" checkbox');
end;
disp('Checking completed.');

author: PatternRecognition

repo: OpenBMI


text:
function orthogonalslicer(imaVOL, pixsize, colormaptype, minpix, maxpix)
% ORTHOGONALSLICER
%
% Orthogonalslicer creates 3 figures showing the 3 orthogonal slices
% from a scalar volume of data. Cross lines on each image represents lines of
% intersection of the appropriate plane and the two other orthogonal plane.
% By left dragging the mouse on any plane the user can explore the volume.
% Additional mouse options:
%   - the images can be zoomed and paned with middle and right mouse dragging
%   - to restore the original state (zoom and pan) double click on the image
%
% Usage:
%
%   orthogonalslicer; or orthogonalslicer(1);
%   DEMO using the MATLAB 'MRI' data set.
%
%   orthogonalslicer(2);
%   DEMO using the MATLAB 'FLOW' data set.
%
%   orthogonalslicer(VOLUME)
%       - VOLUME a 3D scalar data array
%
%   orthogonalslicer(VOLUME, PIXELSIZE)
%       - PIXELSIZE describes the size(aspect) of the pixel in 3D /default = [1 1 1]/
%
%   orthogonalslicer(VOLUME, PIXELSIZE, COLORMAPTYPE)
%       - COLORMAPTYPE selects the colormap to be used /default = jet/
%
%   orthogonalslicer(VOLUME, PIXELSIZE, COLORMAPTYPE,MINPIX,MAXPIX))
%       - MINPIX,MAXPIX normalizes the values in VOLUME to the range
%       specified (same as IMAGESC does)
%

%
% University of Debrecen, PET Center/LB 2004

if nargin == 0
    mridata = load('mri');
    imaVOL = squeeze(mridata.D);
    minpix = 0; maxpix = max(imaVOL(:));
    pixsize = [1 1 1/0.4];
    colormaptype = bone(256);
elseif nargin == 1
    if sum(size(imaVOL)) == 2
        if imaVOL == 1
            mridata = load('mri');
            imaVOL = squeeze(mridata.D);
            minpix = 0; maxpix = max(imaVOL(:));
            pixsize = [1 1 1/0.4];
            colormaptype = bone(256);
        elseif imaVOL == 2
            imaVOL = permute(flow(50),[1 3 2]);
            pixsize = [1 1 1];
            colormaptype = 'jet';
            minpix = min(imaVOL(:)); maxpix = max(imaVOL(:));
        end
        mia_Start3dCursor(imaVOL, pixsize, colormaptype, minpix, maxpix);
        return;
    end
    pixsize = [1 1 1];
    colormaptype = jet;
    minpix = min(imaVOL(:)); maxpix = max(imaVOL(:));
elseif nargin == 2
    colormaptype = jet;
    minpix = min(imaVOL(:)); maxpix = max(imaVOL(:));
elseif nargin == 3
    minpix = min(imaVOL(:)); maxpix = max(imaVOL(:));
elseif nargin == 4
    maxpix = max(imaVOL(:));
end

mia_Start3dCursor(imaVOL, pixsize, colormaptype, minpix, maxpix);

author: Sable

repo: mcbench-benchmarks


text:
Fs=44100; % sampling frequency
dt=1/Fs;
T0=2; % period, sec

close all;

r = audiorecorder(Fs, 16, 1);
clc;disp('recording started...');
recordblocking(r,T0); % record next data
disp('recorded');
s00 = getaudiodata(r); % get data

L0=length(s00);

s=s00-mean(s00);

 %load('bp.mat'); % bandpass 50-4000 Hz
 load('lp.mat'); % lowpass 0-10000Hz

 s = filter(Num,1,s);

 s=s/max(abs(s)); % normalize to maximum

[op ismax]=find_all_optimums(s);
% op - optimums positions

%[z isrise]=find_zeros(s);
%op=decimate_optimums(z,op,s(op));

t=(0:L0-1)'*dt;
plot(t,s,'b-');
hold on;
plot(t(op),s(op),'rx');

% sw0=interp1(t(op),s(op),t,'nearest');
% sw1=interp1(t(op),s(op),t,'linear');
% al=0;

dop=diff(op);

% quntitize dop to [0 255]:
dop(dop>255)=255;
dop=round(dop); % not necessary

opv=s(op); % values in optimums

% quntitize opv to -128 127
opv=128*opv;
opv(opv>127)=127;
opv(opv<-128)=-128;
opv=round(opv);

% compression ratio:
szo=2*length(s); % original size as 2 byts per sample
szc=1*length(dop)+1*length(opv); % compressed size one byte for dop element, one byte for opv element
cr=szo/szc; % compression ratio

% reconstruction
opr=[0; cumsum(dop)];
tr=opr*dt;

opvr=opv/128;

tra=0:dt:max(tr); % time for interpolation

sw=interp1(tr,opvr,t,'pchip');

soundsc(sw,Fs);

plot(t,sw,'g-');

title(['compressed in ' num2str(cr) ' times']);
legend('original','optimums','reconstructed');


text:
%% Efficacy of coadministration of Drug X with Statin on cholesterol reduction
% Copyright 2010 - 2011 MathWorks, Inc.

%% Abstract
% Statins are the most common class of drugs used for treating
% hyperlipdemia. However, studies have shown that even at their maximum
% dosage of 80 mg, many patients do not reach LDL cholesterol goals
% recommended by the National Cholesterol Education Program Adult Treatment
% Panel. Combination therapy, in which a second cholesterol-reducing agent that
% acts via a complementary pathway is coadmininstered with statin, is one
% alternative of achieving higher efficacy at lower statin dosage.
%
% In this example, we test the primary hypothesis that coadminstering drug
% X with statin is more effective at reducing cholesterol levels than
% statin monotherapy.
%
% *NOTE The dataset used in this example is purely fictitious*.
%
% The analysis presented in this example is adapted from the following
% publication.
%
% *Reference* Ballantyne CM, Houri J, Notarbartolo A, Melani L, Lipka LJ,
% Suresh R, Sun S, LeBeaut AP, Sager PT, Veltri EP; _Ezetimibe Study Group.
% Effect of ezetimibe coadministered with atorvastatin in 628 patients with
% primary hypercholesterolemia: a prospective, randomized, double-blind
% trial._ Circulation. 2003 May 20;107(19):2409-15.

%% Data
% 650 patients were randomly assigned to one of the following 10 treatment
% groups (65 subjects per group)
%
% * Placebo
% * Drug X (10 mg)
% * Statin (10, 20, 40 or 80 mg)
% * Drug X (10 mg) + Statin (10, 20, 40 or 80 mg)
%
% Lipid profile (LDL cholesterol, HDL CHolesterol and Triglycerides) was
% measured at baseline (BL) and at 12 weeks (after the start of treatment).
% In addition to the lipid profile, patients age, gender and Cardiac Heart
% Disease (CHD) risk category was also logged at baseline.
%
% The data from the study is stored in a Microsoft Excel (R) file. Note
% that the  data could also be imported from other sources such as text
% files, any JDBC/ODBC compliant database, SAS transport files, etc.
%
% The columns in the data are as follows:
%
% * ID     - Patient ID
% * Group  - Treatment group
% * Dose_A - Dosage of Statin (mg)
% * Dose_X - Dosage of Drug X (mg)
% * Age    - Patient Age
% * Gender - Patient Gender
% * Risk   - Patient CHD risk category (1 is high risk, and 3 is low risk)
% * LDL_BL - HDL_BL & TC_BL - Lipid levels at baseline
% * LDL_12wks , HDL_12wks & TC_12wks - Lipid levels after treatment

%%
% We will import the data into a dataset array that affords better data
% managemment and organization.

% Import data from an Excel file
  ds = dataset('xlsfile', 'Data.xls') ;

%% Preliminary analysis
% Our primary efficacy endpoint is the level of LDL cholesterol. Let us
% compare the LDL C levels at baseline to LDL C levels after treatment

% Use custom scatter plot
  LDLplot(ds.LDL_BL, ds.LDL_12wk, 50, 'g')

%%
% The mean LDL C level at baseline is around 4.2 and mean level after
% treatment is 2.5. So, at least for the data pooled across all the
% treatment groups, it seems that the treatment causes lowering of the LDL
% cholesterol levels

% Use a grouped scatter plot
  figure
  gscatter(ds.LDL_BL, ds.LDL_12wk, ds.Group)

%%
% The grouped plot shows that LDL C levels before the start of treatment
% have similar means. However, the LDL C levels after treatment show
% difference across treatment groups. The Placebo group show no
% improvement. Statin monotherapy seems to outperform the Drug X
% monotherapy. There is overlap between the Statin and Statin + X groups;
% however, it the combination treatment does seem to perform better that
% the statin monotherapy. Remember that the "Statin" and "Statin + X"
% groups are further split based on Statin dose.
%
% In this example, we will use percentage change of LDL C from the baseline
% level as the primary metric of efficacy.

% Calculate the percentage improvement over baseline level

  ds.Change_LDL = ( ds.LDL_BL - ds.LDL_12wk ) ./ ds.LDL_BL * 100 ;

%%
% In the following graph, we can see that
%
% # In the "Statin" and "Statin + X" group, there appears to be a positive
% linear correlation between percentage improvement and statin dose
% # Even at the smallest dose of 10 mg, monotherapy with statin seems to be
% better than the Drug X monotherapy group

% Visualize effect of treatment and statin dose on perecentage LDL reduction
  figure
  gscatter(ds.ID, ds.Change_LDL, {ds.Group, ds.Dose_S})
  legend('Location', 'Best')

%% Pooled comparison: Is the combination therapy better than statin monotherapy ?
%
% First, we will extract percent change in LDL C level for the Statin and
% the Statin + X groups only.  We will test the null hypothesis that the
% percent change in LDL C level for the "Statin + X" groups is greater than
% that in the "Statin + X" using pooled data. We use a 2 sample t-test to
% test this hypothesis.

% Convert Group into a categorical variable
ds.Group = nominal(ds.Group) ;

grp1 = ds.Change_LDL(ds.Group == 'Statin')         ;
grp2 = ds.Change_LDL(ds.Group == 'Statin + X')     ;

[h, p] = ttest2(grp1, grp2, .01, 'left')

%%
% We performed a tailed hypothesis to see if Statin + X group (grp2) is
% better than the Statin group (grp1). We test against the alternative that
% that mean LDL change of grp1 (Statin only) is less than mean LDL change
% of grp2 (Statin + X)
%
% The null hypothesis is rejected (p < 0.01), implying that grp1 mean is
% less that grp2 mean, i.e. the Statin group is less effective at lowering
% LDL C levels than the Statin + X group.
%
% The pooled analysis shows that coadministering drug X with statin is more
% effective than statin monotherapy.

%% Effect of Treatment, Statin Dose and Dose by Treatment interaction
% Our analysis so far was done on pooled data. We analysed the effect of
% treatment (statin alone (X = 0) vs. statin + 10 mg X) on the LDL C
% levels. We ignored levels of statin dose within each treatment group
%
% Next, we will perform a 2-way ANOVA (analysis of variance) to
% simultaneously understand the effect of both factors - statin dose (4
% levels -  10 20, 40, 80 mg) and Treatment (2 level - statin only or
% Statin + 10 mg X ) - on the percentage change of LDL C levels.
%

% First, we filter the data to include only the Statin and Statin + X groups
ds1       = ds(ds.Group == 'Statin' | ds.Group == 'Statin + X', :)    ;

anovan(ds1.Change_LDL , {ds1.Dose_S, ds1.Group }         , ...
        'varnames'    , {'Statin Dose', 'Treatment'}    ) ;

%% Effect of Statin Dose on incremental increase in percentage LDL reduction
% The ANOVA results indicate that statin dose is a significant factor, but
% it doesn't compare means across individual dose-treatment level
% combination. Let's look at the individual cell means.

ds2 = grpstats(ds1 , {'Dose_X', 'Dose_S'}, '', 'DataVars', 'Change_LDL')

%%
% Convert to wide format
ds2 = unstack(ds2, 'mean_Change_LDL' , 'Dose_X', ...
                   'NewDataVarNames' , {'Change_LDL_St', 'Change_LDL_St_X'} )

%%
% From the above table, we can clearly see that the average efficacy of the
% combination therapy is better than statin monotherapy at all statin
% dosages.
%
% In the plot of the individual means, notice that the percentage reduction
% in LDL C levels achieved in the low dose combination therapy group (~50.5
% %) is comparable to that achieved in the higher dose Statin monotherapy
% group (~ 49.4 %). Thus combination therapy with Drug X could help
% patients that cannot tolerate high statin doses.

figure
bar([ds2.Change_LDL_St, ds2.Change_LDL_St_X])
set(gca, 'XTickLabel', [10, 20 40, 80])
colormap summer
xlabel('Statin Dose Groups(mg)')
ylabel('Percentage reduction of LDL C from Baseline (mmol/L)')
legend('Statin', 'Statin + X')

%% Regression analysis: Effect of statin dose on percent LDL C reduction
% In the above graph, there appears to be linear improvement in the
% effectiveness metric for both treatment groups. In general it seems that
% for every doubling of the statin dose, there is a 5-6 point improvement
% in the percentage LDL C reduction. Let's fit a linear regression line to
% the entire dataset, instead of to the mean level.
%
x = ds1.Dose_S    (  ds1.Group == 'Statin' ) ;
y = ds1.Change_LDL(  ds1.Group == 'Statin' ) ;

x1 = ds1.Dose_S    (ds1.Group == 'Statin + X')  ;
y1 = ds1.Change_LDL(ds1.Group == 'Statin + X')  ;

%%

cftool

%%
% The regression line for the Statin and the Statin + X group run almost
% parallel. This probably indicates mechanism of actions of drug X and
% statins are independent.

% Fit
[m1, m2] = createFit(x,y,x1, y1)

%% Secondary Analysis: Consistency of effect across subgroups, age and gender
% Finally, we will make a visual check to ensure that the efficacy of the
% Statin + X treatment at various statin doses is consistent across gender
% and age subgroups. We will perform this check for only the Statin + X
% treatment group.

idx = ds.Group == 'Statin + X' ;
boxplot(ds.Change_LDL(idx), { ds.Dose_S(idx), ds.Gender(idx)} )

%%
% We will convert the continuous age variable into a catergorical variable,
% with 2 categories: Age < 65 and Age >= 65

% Convert age into a ordinal array
ds.AgeGroup = ordinal(ds.Age ,{'< 65', '>= 65'} , [] ,[0 65 100] )  ;

% Plot
boxplot(ds.Change_LDL(idx), { ds.Dose_S(idx), ds.AgeGroup(idx)} )


text:
function out=acosh(x)

out=log(x+sqrt(x+1).*sqrt(x-1));

%%%precision=x(1).precision;
%%%out=mp(zeros(size(x)));
%%%
%%%mpPi=mppi(precision);
%%%for ii=1:numel(x)
%%% out(ii)=log(x(ii)+sqrt(x(ii)+1)*sqrt(x(ii)-1));
%%%end % for ii=1:max(ex,


text:
function prob_test1188 ( )

%*****************************************************************************80
%
%% PROB_TEST1188 tests NORMAL_TRUNCATED_B_CDF, NORMAL_TRUNCATED_B_CDF_INV, NORMAL_TRUNCATED_B_PDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 August 2013
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;
  b = 150.0;
  mu = 100.0;
  s = 25.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROB_TEST1188\n' );
  fprintf ( 1, '  For the Upper Truncated Normal PDF:\n' );
  fprintf ( 1, '  NORMAL_TRUNCATED_B_CDF evaluates the CDF.\n' );
  fprintf ( 1, '  NORMAL_TRUNCATED_B_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  NORMAL_TRUNCATED_B_PDF evaluates the PDF.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The "parent" normal distribution has\n' );
  fprintf ( 1, '    mean =               %g\n', mu );
  fprintf ( 1, '    standard deviation = %g\n', s );
  fprintf ( 1, '  The parent distribution is truncated to\n' );
  fprintf ( 1, '  the interval [-oo,%g]\n', b );

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

  for i = 1 : 10

    [ x, seed ] = normal_truncated_b_sample ( mu, s, b, seed );

    pdf = normal_truncated_b_pdf ( x, mu, s, b );

    cdf = normal_truncated_b_cdf ( x, mu, s, b );

    x2 = normal_truncated_b_cdf_inv ( cdf, mu, s, b );

    fprintf( 1, '  %14g  %14g  %14g  %14g\n', x, pdf, cdf, x2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear;
data = load('data_completeness.mat');
data = data.data;

orange = [ 0.91,0.41,0.17];
blue = [0,0.4470,0.7410];
figure('DefaultAxesFontSize',14);
ax = gca;
xlabel('Month');

yyaxis left;
h1 = plot(data(:,1),'-','color',blue,'linewidth',2,'DisplayName','rmse');
ylabel('Data amount (\times 10^6)');
ylim([0,1.6e+6]);
ax.YTick = 0:0.4e+6:1.6e+6;
ax.YTickLabel = ['0  ';'0.4';'0.8';'1.2';'1.6'];
yyaxis right;
h2 = plot(data(:,2),'-','color',orange,'linewidth',2,'DisplayName','lb');
ylabel('Completeness');
ax.YTick = 0:0.25:1;
ax.YTickLabel = 0:0.25:1;

box on;
grid on;
ax.XTick = 1:4:48;
ax.XTickLabel = ['Apr.';'Aug.';'Dec.';'Apr.';'Aug.';'Dec.';'Apr.';'Aug.';'Dec.';'Apr.';'Aug.';'Dec.';];
ax.GridLineStyle = '-.';
set(gcf, 'PaperSize', [8 3]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 8 3]);
saveas(gcf,'nyc_data_completeness','pdf');

author: xinychen

repo: academic-drawing


text:
function knapsack_01_test ( )

%*****************************************************************************80
%
%% KNAPSACK_01_TEST tests the KNAPSACK_01 library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 August 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'KNAPSACK_01_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the KNAPSACK_01 library.\n' );

  knapsack_01_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'KNAPSACK_01_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% Copyright 2017 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     https://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

% Returns set of parameters
% set light_mode = true to run the code in a sub optimal but faster mode
% set light_mode = false to obtain the results reported in the RED paper

function params = GetSuperResSDParams(light_mode)

% regularization factor
params.lambda = 0.008;

% number of outer iterations
if light_mode
    params.outer_iters = 500;
else
    params.outer_iters = 1500;
end

% level of noise assumed in the regularization-denoiser
params.effective_sigma = 3;

return

author: google

repo: RED

save_path: github-repos/MATLAB/google-RED