text:
function pce_burgers ( )

%*****************************************************************************80
%
%% PCE_BURGERS applies the polynomial chaos expansion to the Burgers equation.
%
%  Discussion:
%
%    The time-dependent viscous Burgers equation to be solved is:
%
%      du/dt = - d ( u*(1/2-u)) /dx + nu d2u/dx2
%
%    with boundary conditions
%
%      u(-3.0) = 0.0, u(+3.0) = 1.0.
%
%    The viscosity nu is assumed to be an uncertain quantity with
%    normal distribution of known mean and variance.
%
%    A polynomial chaos expansion is to be used, with Hermite polynomial
%    basis functions h(i,x), 0 <= i <= n.
%
%    Because the first two Hermite polynomials are simply 1 and x,
%    we have that
%
%      nu = nu_mean * h(0,x) + nu_variance * h(1,x).
%
%    We replace the time derivative by an explicit Euler approximation,
%    so that the equation now describes the value of U(x,t+dt) in terms
%    of known data at time t.
%
%    Now assume that the solution U(x,t) can be approximated
%    by the truncated expansion:
%
%      U(x,t) = sum ( 0 <= i <= n ) c(i,t) * h(i,x)
%
%    In the equation, we replace U by its expansion, and then multiply
%    successively by each of the basis functions h(*,x) to get a set of
%    n+1 equations that can be used to determine the values of c(i,t+dt).
%
%    This process is repeated until the desired final time is reached.
%
%    At any time, the coefficients c(0,t) contain information definining
%    the expected value of u(x,t) at that time, while the higher order coefficients
%    can be used to deterimine higher moments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2012
%
%  Author:
%
%    Original FORTRAN90 version by Gianluca Iaccarino.
%    This MATLAB version is by John Burkardt.
%
%  Local parameters:
%
%    Local, real DT, the timestep.
%
%    Local, real DX, the spacing between grid points.
%
%    Local, integer N, the number of intervals in the spatial domain.
%
%    Local, real NUMEAN, the mean of viscosity.
%
%    Local, real NUVARIANCE, the variance of viscosity.
%
%    Local, integer P, the order of the PC expansion.
%
%    Local, real T, the current time.
%
%    Local, real TF, the final integration time.
%
%    Local, real U1(N+1,P+1), the PCE representation at the current time.
%
%    Local, real U2(N+1,P+1), the PCE representation for the next time.
%
%    Local, real X(N+1,1), the grid points.
%
  p = 5;
  n = 32;
  nt = 2000;
  ti = 0.0;
  tf = 2.0;
  dt = ( tf - ti ) / nt;
  numean = 0.25;
  nuvariance = 0.08;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PCE_BURGERS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Polynomial Chaos Expansion\n' );
  fprintf ( 1, '  1D Burgers equation\n' );
  fprintf ( 1, '  Original version by Gianluca Iaccarino\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  PCE order       = %d\n', p );
  fprintf ( 1, '  Number of cells = %d\n', n );
  fprintf ( 1, '  Time step       = %g\n', dt );
  fprintf ( 1, '  Initial time    = %g\n', ti );
  fprintf ( 1, '  Final time      = %g\n', tf );
  fprintf ( 1, '  Viscosity Mean  = %g\n', numean );
  fprintf ( 1, '  Viscosity Var   = %g\n', nuvariance );
  fprintf ( 1, '\n' );
%
%  Define some numerical parameters.
%
  dx = 6.0 / n;
  conv = dt / ( 2.0 * dx );

  visc = zeros ( p + 1, 1 );
  visc(1) = numean * dt / ( dx * dx );
  visc(2) = nuvariance * dt / ( dx * dx );
%
%  Define a uniform grid.
%
  x = ( linspace ( - 3.0, + 3.0, n + 1 ) )';
%
%  Set the initial conditions.
%
  u1 = zeros ( n + 1, p + 1 );
  u1(1:n+1,1) = 0.5 + x(1:n+1,1) / 6.0;

  u2 = zeros ( n + 1, p + 1 );
%
%  Time integration.
%
  t1 = ti;
%
%  Write the current solution.
%
  output_filename = 'burgers.history.txt';
  output_unit = fopen ( output_filename, 'wt' );
  fprintf ( output_unit, '----------\n' );
  fprintf ( output_unit, 'T = %g\n', t1 );
  for i = 1 : n + 1
    fprintf ( output_unit, '  %10g', x(i) );
    for k = 1 : p + 1
      fprintf ( output_unit, '  %10g', u1(i,k) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Time integration.
%
  for it = 1 : nt

    t2 = ( ( nt - it ) * ti   ...
         + (      it ) * tf ) ...
         / ( nt      );
%
%  Boundary conditions.
%
    u2(1,1:p+1) = 0.0;
    u2(n+1,1) = 1.0;
    u2(n+1,2:p+1) = 0.0;

    for k = 1 : p + 1

      dp = he_double_product_integral ( k - 1, k - 1 );

      for ix = 2 : n
%
%  Viscous term.
%
        term1 = visc(1) * ( u1(ix+1,k) - 2.0 * u1(ix,k) + u1(ix-1,k) );
        i = 2;
        for j = 1 : p + 1
          tp = he_triple_product_integral ( i - 1, j - 1, k - 1 );
          term1 = term1 + visc(i) * ( u1(ix+1,j) - 2.0 * u1(ix,j) + u1(ix-1,j) ) * tp / dp;
        end
%
%  Convective term.
%
        term2 = - conv * 0.5 * ( u1(ix+1,k) - u1(ix-1,k) );
        for j = 1 : p + 1
          for i = 1 : p + 1
            tp = he_triple_product_integral ( i - 1, j - 1, k - 1 );
            term2 = term2 + ( conv * u1(ix,i) * ( u1(ix+1,j) - u1(ix-1,j) ) * tp ) / dp;
          end
        end

        u2(ix,k) = u1(ix,k) + term1 + term2;

      end

    end

    t1 = t2;
    u1(1:n+1,1:p+1) = u2(1:n+1,1:p+1);
%
%  Print solution every 100 time steps.
%
    if ( mod ( it, 100 ) == 0 )
      fprintf ( output_unit, '----------\n' );
      fprintf ( output_unit, 'T = %g\n', t1 );
      for i = 1 : n + 1
        fprintf ( output_unit, '  %10g', x(i) );
        for k = 1 : p + 1
          fprintf ( output_unit, '  %10g', u1(i,k) );
        end
        fprintf ( output_unit, '\n' );
      end
    end

  end

  fclose ( output_unit );
  fprintf ( 1, '  Time history in "%s".\n', output_filename );
%
%  Compute the mean and variance.
%
  umean(1:n+1,1) = u1(1:n+1,1);

  uvariance = zeros ( n + 1, 1 );

  for i = 1 : n + 1
    for j = 2 : p + 1
      dp = he_double_product_integral ( j - 1, j - 1 );
      uvariance(i) = uvariance(i) + u1(i,j).^2 * dp;
    end
  end
%
%  Save the solution at the final time.
%
  output_filename = 'burgers.moments.txt';
  output_unit = fopen ( output_filename, 'wt' );
  fprintf ( output_unit, ' X E[U] Var[U]\n' );
  for i = 1 : n + 1
    fprintf ( output_unit, '  %18.8g  %18.8g  %18.8g\n', ...
      x(i), umean(i), uvariance(i) );
  end
  fclose ( output_unit );
  fprintf ( 1, '  Moments in "%s".\n', output_filename );

  output_filename = 'burgers.modes.txt';
  output_unit = fopen ( output_filename, 'wt' );
  fprintf ( output_unit, ' X U_0 ... U_P \n' );
  for i = 1 : n + 1
    fprintf ( output_unit, '  %10g', x(i) );
    for k = 1 : p + 1
      fprintf ( output_unit, '  %10g', u1(i,k) );
    end
    fprintf ( output_unit, '\n' );
  end
  fclose ( output_unit );
  fprintf ( 1, '  Final modes in "%s".\n', output_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PCE_BURGERS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function jed = epoch_to_jed_ethiopian ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_ETHIOPIAN returns the epoch of the Ethiopian calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 1724220.5;

  return
end


text:
function K = rbf(kern,dat1,dat2,ind1,ind2,kerParam),

% K = rbf(d1,d2,ind1,ind2,param), compute the kernel
%     matrix between d1 and d2
% for a rbf kernel exp(-||x-z||^2/(2*param^2))
%     where x is from d1 and z from d2

  K=get_x(dat2,ind2)*get_x(dat1,ind1)';
  kernTemp=kernel;
  Kdn = get_norm(kernTemp,dat1,ind1).^2;
  Kn = get_norm(kernTemp,dat2,ind2).^2;
  K = ones(length(Kn),1)*Kdn' + Kn*ones(1,length(Kdn)) - 2*K;
  K = exp(-K/(2*kerParam^2));

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef Triangle_Linear < Interpolation

    methods (Access = public)

        function obj = Triangle_Linear(cParams)
            obj.init(cParams);
            obj.computeParams();
            obj.computeCases();
        end

        function computeShapeDeriv(obj,xGauss)
            obj.computeShapes(xGauss);
            obj.computeShapeDerivatives(xGauss);
        end

    end

    methods (Access = private)

        function computeParams(obj)
            obj.ndime = 2;
            obj.nnode = 3;
            obj.pos_nodes = [0 0; 1 0; 0 1];
            obj.isoDv = 0.5;
            obj.main_loop = [3 3];
            obj.extra_cases = [];
        end

        function computeShapes(obj,posgp)
            ngaus = size(posgp,2);
            nelem = size(posgp,3);
            s = posgp(1,:,:);
            t = posgp(2,:,:);
            I = ones(size(t));
            obj.shape = zeros(obj.nnode,ngaus,nelem);
            obj.shape(1,:,:) = I-s-t;
            obj.shape(2,:,:) = s;
            obj.shape(3,:,:) = t;
        end

        function computeShapeDerivatives(obj,posgp)
            ngaus = size(posgp,2);
            nelem = size(posgp,3);
            obj.deriv = zeros(obj.ndime,obj.nnode,ngaus,nelem);
            obj.deriv(1,1,:,:) = -1;
            obj.deriv(1,2,:,:) = 1;
            obj.deriv(1,3,:,:) = 0;
            obj.deriv(2,1,:,:) = -1;
            obj.deriv(2,2,:,:) = 0;
            obj.deriv(2,3,:,:) = 1;
        end

        function computeCases(obj)
            obj.iteration = [1 2 3;
                             2 3 1];
            obj.cases(:,:,1) = [1 4 5;
                4 2 3;
                5 4 3];
            obj.cases(:,:,2) = [1 4 3;
                4 2 5;
                4 5 3];
            obj.cases(:,:,3) = [1 4 5;
                1 2 4;
                5 4 3];
            obj.selectcases =  [1     0;
                2     0;
                3     3;
                0     2;
                0     1];
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function rbfBasis = rbfComputeBasis( X, k, cluster, scale, show )
% Get locations and sizes of radial basis functions for use in rbf network.
%
% Radial basis function are a simple, fast method for universal function
% approximation / regression. The idea is to find a mapping X->y where X
% and y are continuous real variables. The mapping is linear over features
% of X: y=rbfWeight*features(X). The features are of the form:
%  f_j(x) = exp(-||x-mu_j||_2^2 / 2 sig_j^2 ).
% The number of basis functions controls the complexity of the network. Too
% many basis functions can lead to overfitting, too few to bad
% interpolation. Rbf networks are trained in two phases.  First, the radial
% basis functions are found in an unsupervised manner (using only Xtrain),
% for example by clustering Xtrain.  Next, given the basis functions,
% Xtrain and ytrain, the basis weights are found by solving the system:
%  rbfWeight * features(Xtrain) = ytrain
% At this point, to interpolate any new points, Xtest, use:
%  ytest = rbfWeight * features(Xtest)
% The code below achieves all three steps:
%  rbfBasis  = rbfComputeBasis( Xtrain, k, cluster, scale, show );
%  rbfWeight = rbfComputeFtrs(Xtrain,rbfBasis) \ ytrain;
%  ytest     = rbfComputeFtrs(Xtest,rbfBasis) * rbfWeight;
% Note, in the returned rbfBasis struct there are a number of flags that
% control how the rbf features are computed. These can be altered to
% achieve the desired effect.
%
% For an in depth discussion of rbf networks see:
%  Christopher M. Bishop. "Neural Networks for Pattern Recognition"
%
% USAGE
%  rbfBasis = rbfComputeBasis( X, k, [cluster], [scale], [show] )
%
% INPUTS
%  X           - [N x d] N points of d dimensions each
%  k           - number of basis functions to use
%  cluster     - [1]: Computes cluster centers for use as rbf functions.
%              - 0: Evenly centered basis functions (ok for small d)
%  scale       - [5] Alter computed value of sigma by given factor
%                set larger for smoother results, too small -> bad interp
%  show        - [0] will display results in figure(show)
%                if show<0, assumes X is array Nxs^2 of N sxs patches
%
% OUTPUTS
%  rfbBasis
%   .d          - feature vector size
%   .k          - number of basis functions actually used
%   .mu         - [d x k] rbf centers
%   .vars       - [1 x k] rbf widths
%   .var        - rbf average width
%   .globalVar  - [1] if true use single average var for rbfs
%   .constant   - [0] if true include extra basis with constant activation
%   .normalize  - [0] if true normalize overall rbf response to sum to 1
%
% EXAMPLE
%
% See also RBFDEMO, RBFCOMPUTEFTRS
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.50
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<2 || isempty(k)); error('k not specified'); end
if( nargin<3 || isempty(cluster)); cluster=1;  end
if( nargin<4 || isempty(scale)); scale=5;  end
if( nargin<5 || isempty(show)); show=0;  end
[N, d] = size(X);

if( cluster )
  %%% CLUSTERS subsample, run kmeans
  maxN=5000; if( N>maxN );  X=X(randSample(N,maxN),:);  N=maxN;  end
  prm.nTrial=5; prm.display=0;
  [IDX,mu] = kmeans2( X, k, prm );
  mu = mu'; k = size(mu,2);
else
  %%% GRID generate locations evenly spaced on grid
  if( d>4 ); error('d too high. curse of dimensionality..'); end
  nBPer = round( k ^ (1/d) );  k = nBPer ^ d; rg=[min(X)' max(X)'];
  del=(rg(:,2)-rg(:,1))/(nBPer-1); rg=rg+[-del del]/2;
  loc=cell(1,d);  for i=1:d; loc{i}=linspace(rg(i,1),rg(i,2),nBPer); end
  grid=cell(1,d); if(d>1); [grid{:}]=ndgrid(loc{:}); else grid=loc; end
  mu=zeros(d,k);   for i=1:d; mu(i,:) = grid{i}(:); end
end

%%% Set var to be equal to average distance of neareast neighbor.
dist = pdist2( mu', mu' );
dist = dist + realmax * eye( k );
vars = min(dist)* scale;
var  = mean(vars);
vars = max( vars, var/100 );

%%% store results
rbfBasis = struct('d',d, 'k',k, 'mu',mu, 'vars',vars, 'var',var, ...
  'globalVar',1, 'constant',0, 'normalize',0);

%%% optionally display
if( abs(show) )
  if( show<0 ) % if images can display
    siz = sqrt(d);
    I = clusterMontage( reshape(X,siz,siz,N), IDX, 25, 1 );
    figure(-show); clf; montage2( I );
    figure(-show+1); clf; montage2(reshape(mu,siz,siz,[]));
  elseif( d==1 ) % 1D data
    figure(show); clf; hold on;
    minX = min(X,[],1 );  maxX = max(X,[],1 );
    xs = linspace( minX, maxX, 500 )';
    for i=1:k
      ys = exp( -(xs-mu(i)).^2 / 2 / var  );
      plot( xs, ys );
    end
  elseif( d==2 ) % 2D data
    figure(show); clf;
    minX = min(X,[],1 );  maxX = max(X,[],1 );
    xs1 = linspace(minX(1),maxX(1),25);
    xs2 = linspace(minX(2),maxX(2),25);
    [xs1,xs2] = ndgrid( xs1, xs2 );
    xs = [xs1(:) xs2(:)]; n = size(xs,1);
    for i=1:k
      mui = repmat(mu(:,i),[1 n])';
      ys = exp( - sum( ((xs - mui)).^2, 2 ) / 2 / var );
      surf( xs1, xs2, reshape(ys,size(xs1)) );
      hold on;
    end;
  elseif( d==3 ) % 3D data (show data+centers)
    figure(show); clf; hold on;
    scatter3( X(:,1),X(:,2),X(:,3),12,'filled');
    scatter3( mu(1,:),mu(2,:),mu(3,:),30,'filled');
  end
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function p=protate(p,phi)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

A=[cos(phi),-sin(phi);sin(phi),cos(phi)];
p=p*A;

repo: ifem


text:
function [phi, lam, h, phiC] = cart2geod(X, Y, Z)

% SYNTAX:
%   [phi, lam, h, phiC] = cart2geod(X, Y, Z);
%   [phi, lam, h, phiC] = cart2geod(xyz);
%
% INPUT:
%   X = X axis cartesian coordinate
%   Y = Y axis cartesian coordinate
%   Z = Z axis cartesian coordinate
%
% OUTPUT:
%   phi = latitude
%   lam = longitude
%   h = ellipsoidal height
%   phiC = geocentric latitude
%
% DESCRIPTION:
%   Conversion from cartesian coordinates to geodetic coordinates.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%global a_GPS e_GPS

if nargin == 1
    Z = X(:, 3);
    Y = X(:, 2);
    X = X(:, 1);
end

a = GPS_SS.ELL_A;
e = GPS_SS.ELL_E;

%radius computation
r = sqrt(X.^2 + Y.^2 + Z.^2);

%longitude
lam = atan2(Y,X);

%geocentric latitude
phiC = atan(Z./sqrt(X.^2 + Y.^2));

%coordinate transformation
psi = atan(tan(phiC)/sqrt(1-e^2));

phi = atan((r.*sin(phiC) + e^2*a/sqrt(1-e^2) * (sin(psi)).^3) ./ ...
    			(r.*cos(phiC) - e^2*a * (cos(psi)).^3));

N = a ./ sqrt(1 - e^2 * sin(phi).^2);

%height
h = r .* cos(phiC)./cos(phi) - N;

author: goGPS-Project


text:
classdef UnfittedMeshPlotter < handle

    properties (Access = private)
        uMesh
        faceColor
        edgeAlpha
        faceAlpha
    end

    methods (Access = public)

        function obj = UnfittedMeshPlotter(cParams)
            obj.init(cParams)
        end

        function plotDomainInColor(obj,color)
            obj.faceColor = color;
            obj.edgeAlpha = 0;
            obj.faceAlpha = 1;
            obj.plotInner();
            obj.plotInnerCut();
        end

        function plotDomain(obj)
            obj.faceColor = [];
            obj.edgeAlpha = [];
            obj.faceAlpha = [];
            obj.plotAll();
            obj.addLighting();
        end

        function plotBoundary(obj)
            obj.plotBackground();
            obj.plotBoundaryCutMesh();
            obj.plotUnfittedBoundaryMesh();
            obj.addLighting();
        end

        function plotAll(obj)
            obj.plotBackground();
            obj.plotInner();
            obj.plotInnerCut();
            obj.plotBoundaryCutMesh();
            obj.plotUnfittedBoundaryMesh();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.uMesh = cParams.uMesh;
        end

        function plotBackground(obj)
            s.mesh = obj.uMesh.backgroundMesh;
            s.isBackground = true;
            obj.plotMesh(s);
        end

        function plotInner(obj)
            m = obj.uMesh.innerMesh;
            obj.plotSubMesh(m);
        end

        function plotInnerCut(obj)
            m = obj.uMesh.innerCutMesh;
            obj.plotSubMesh(m);
        end

        function plotBoundaryCutMesh(obj)
            m = obj.uMesh.boundaryCutMesh;
            obj.plotSubMesh(m);
        end

        function plotUnfittedBoundaryMesh(obj)
            uB = obj.uMesh.unfittedBoundaryMesh;
            if ~isempty(uB)
                uBoundaryMeshes = uB.getActiveMesh();
                for imesh = 1:numel(uBoundaryMeshes)
                    uM = uBoundaryMeshes{imesh};
                    uM.plotAll();
                end
            end
        end

        function plotSubMesh(obj,uM)
            if ~isempty(obj.faceColor)
                s.faceColor = obj.faceColor;
            end
            if ~isempty(obj.edgeAlpha)
                s.edgeAlpha = obj.edgeAlpha;
            end
            if ~isempty(obj.faceAlpha)
                s.faceAlpha = obj.faceAlpha;
            end
            if ~isempty(uM)
                s.mesh = uM.mesh;
                obj.plotMesh(s);
            end
        end

        function addLighting(obj)
            delete(findall(gcf,'Type','light'))
            if ~isempty(obj.uMesh.innerMesh)
                if obj.uMesh.innerMesh.mesh.ndim == 3
                    l = lightangle(120,30);
                    material metal
                    l = lightangle(240,30);
                    material metal
                    l = lightangle(90,180);
                    material metal
                    l = lightangle(90,0);
                    material metal
                    view([1 1 1])
                end
            end
        end

    end

    methods (Access = private, Static)

        function plotMesh(s)
            s = SettingsMeshPlotter(s);
            if ~isempty(s.mesh)
                mP = MeshPlotter(s);
                mP.plot();
            end
        end
    end

end


text:
clear;
clc;
format short;

%% 线性代数及其应用第5版 5.8 Iterative Estimates For EienValues

%% EXAMPLE 1
A = [1.8 0.8; 0.2 1.2];
v1 = [4 1]';
x = [-0.5 1]';

for i = 1:8
    x = A*x;
end

x = x / norm(x);
v1= v1 / norm(v1);
%
% fprintf("可以看到 x 逼近v1的方向\r\n");
% x
% v1

%% EXAMPLE 2
A  = [6 5; 1 2];
x = [0 1]';

for i = 1:8
    x = A*x;
    lanbda = abs(max(x));
    x  = x / lanbda;

end

x
lanbda

author: yandld


text:
% Dual RPCA (Lin et al. 2009)
% process_video('RPCA', 'DUAL', 'dataset/demo.avi', 'output/demo_DUAL.avi');
lambda = 1/sqrt(max(size(M))); % default lambda
[L,S] = dual_rpca_2(M,lambda);

author: andrewssobral

repo: lrslibrary


text:
%kPolyline 'Polyline Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros Polyline.pane file
%
% Parameters:
%
% Example: kPolyline( {})
%
% Khoros helpfile follows below:
% .begin tagged
%
% .item "The Polyline Visual Object"
% The polyline object supports the display of a polyline annotation.
% The polyline object is currently still under construction.
%
% .item "Filled"
% When set to true, this attribute indicates that the polyline should be filled
% with the foreground color;  objects are not filled by default.
%
% .item "Border Width"
% This list selection allows you to set the line width used by
% the polyline object. Line widths range from "Extra Fine" to "Extra Wide".
%
% .item "Border Type"
% This list selection lets you set the line type used on
% the polyline object. There are seven line types, including "Solid",
% "Dotted", "Dot Dash", "Short Dash", "Long Dash", "Odd Dash", and
% "Grid Dotted".
%
% .item "Border Color"
% This stringlist selection lets you set the foreground color
% of the polyline object, using the color name.
%
% .item "Fill Color"
% This stringlist selection lets you set the background color
% of the polyline object, using the color name.
% .end tagged

function varargout = kPolyline(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kPolyline(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={};
maxval={};
minval={};
istoggle=[];
was_set=istoggle * 0;
paramtype={};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'Polyline" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [NABF]=analysis_nabf(f,I1,I2)
% function [QABF,LABF,NABF,NABF1]=objective_fusion_perform_fn(f,I1,I2)
%
%%% objective_fusion_perform_fn: Computes the Objective Fusion Performance Parameters proposed by Petrovic
%%% and modified Fusion Artifacts (NABF) measure proposed by B. K. Shreyamsha Kumar
%%%
%%% Inputs:
%%% xrcw -> fused image
%%% x -> source images, x{1}, x{2}
%%%
%%% Outputs:
%%% QABF -> Total information transferred from source images to fused image measure proposed by Petrovic
%%% LABF -> Total loss of information measure proposed by Petrovic
%%% NABF1 -> Fusion Artifacts measure proposed by Petrovic
%%% NABF -> Modified Fusion Artifacts measure proposed by B. K. Shreyamsha Kumar
%%%
%%% Author : B. K. SHREYAMSHA KUMAR
%%% Created on 28-10-2011.
%%% Updated on 08-11-2011.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Petrovic Metrics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Parameters for Petrovic Metrics Computation.
Td=2;
wt_min=0.001;
P=1;
Lg=1.5;
Nrg=0.9999;
kg=19;
sigmag=0.5;
Nra=0.9995;
ka=22;
sigmaa=0.5;

xrcw = (f);
x1 = (I1);
x2 = (I2);

%%% Edge Strength & Orientation.
[gvA,ghA]=sobel_fn(x1);
gA=sqrt(ghA.^2+gvA.^2);

[gvB,ghB]=sobel_fn(x2);
gB=sqrt(ghB.^2+gvB.^2);

[gvF,ghF]=sobel_fn(xrcw);
gF=sqrt(ghF.^2+gvF.^2);

%%% Relative Edge Strength & Orientation.
[p,q]=size(xrcw);
for ii=1:p
   for jj=1:q
      if(gA(ii,jj)==0 | gF(ii,jj)==0)
         gAF(ii,jj)=0;
      elseif(gA(ii,jj)>gF(ii,jj))
         gAF(ii,jj)=gF(ii,jj)/gA(ii,jj);
      else
         gAF(ii,jj)=gA(ii,jj)/gF(ii,jj);
      end
      if(gB(ii,jj)==0 | gF(ii,jj)==0)
         gBF(ii,jj)=0;
      elseif(gB(ii,jj)>gF(ii,jj))
         gBF(ii,jj)=gF(ii,jj)/gB(ii,jj);
      else
         gBF(ii,jj)=gB(ii,jj)/gF(ii,jj);
      end
      if(gvA(ii,jj)==0 & ghA(ii,jj)==0)
         aA(ii,jj)=0;
      else
         aA(ii,jj)=atan(gvA(ii,jj)/ghA(ii,jj));
      end
      if(gvB(ii,jj)==0 & ghB(ii,jj)==0)
         aB(ii,jj)=0;
      else
         aB(ii,jj)=atan(gvB(ii,jj)/ghB(ii,jj));
      end
      if(gvF(ii,jj)==0 & ghF(ii,jj)==0)
         aF(ii,jj)=0;
      else
         aF(ii,jj)=atan(gvF(ii,jj)/ghF(ii,jj));
      end
   end
end
aAF=abs(abs(aA-aF)-pi/2)*2/pi;
aBF=abs(abs(aB-aF)-pi/2)*2/pi;

%%% Edge Preservation Coefficient.
QgAF=Nrg./(1+exp(-kg*(gAF-sigmag)));
QaAF=Nra./(1+exp(-ka*(aAF-sigmaa)));
QAF=sqrt(QgAF.*QaAF);
QgBF=Nrg./(1+exp(-kg*(gBF-sigmag)));
QaBF=Nra./(1+exp(-ka*(aBF-sigmaa)));
QBF=sqrt(QgBF.*QaBF);

%%% Total Fusion Performance (QABF).
wtA=wt_min*ones(p,q);
wtB=wt_min*ones(p,q);
cA=ones(p,q); cB=ones(p,q);
for ii=1:p
   for jj=1:q
      if(gA(ii,jj)>=Td)
         wtA(ii,jj)=cA(ii,jj)*gA(ii,jj)^Lg;
      end
      if(gB(ii,jj)>=Td)
         wtB(ii,jj)=cB(ii,jj)*gB(ii,jj)^Lg;
      end
   end
end
wt_sum=sum(sum(wtA+wtB));
QAF_wtsum=sum(sum(QAF.*wtA))/wt_sum;  %% Information Contributions of A.
QBF_wtsum=sum(sum(QBF.*wtB))/wt_sum;  %% Information Contributions of B.
QABF=QAF_wtsum+QBF_wtsum;   %% QABF=sum(sum(QAF.*wtA+QBF.*wtB))/wt_sum -> Total Fusion Performance.

%%% Fusion Gain (QdeltaABF).
Qdelta=abs(QAF-QBF);
QCinfo=(QAF+QBF-Qdelta)/2;
QdeltaAF=QAF-QCinfo;
QdeltaBF=QBF-QCinfo;
QdeltaAF_wtsum=sum(sum(QdeltaAF.*wtA))/wt_sum;
QdeltaBF_wtsum=sum(sum(QdeltaBF.*wtB))/wt_sum;
QdeltaABF=QdeltaAF_wtsum+QdeltaBF_wtsum;   %% Total Fusion Gain.
QCinfo_wtsum=sum(sum(QCinfo.*(wtA+wtB)))/wt_sum;
QABF11=QdeltaABF+QCinfo_wtsum;              %% Total Fusion Performance.

%%% Fusion Loss (LABF).
rr=zeros(p,q);
for ii=1:p
   for jj=1:q
      if(gF(ii,jj)<=gA(ii,jj) | gF(ii,jj)<=gB(ii,jj))
         rr(ii,jj)=1;
      else
         rr(ii,jj)=0;
      end
   end
end
LABF=sum(sum(rr.*((1-QAF).*wtA+(1-QBF).*wtB)))/wt_sum;

%%% Fusion Artifacts (NABF) by Petrovic.
for ii=1:p
   for jj=1:q
      if(gF(ii,jj)>gA(ii,jj) & gF(ii,jj)>gB(ii,jj))
         na1(ii,jj)=2-QAF(ii,jj)-QBF(ii,jj);
      else
         na1(ii,jj)=0;
      end
   end
end
NABF1=sum(sum(na1.*(wtA+wtB)))/wt_sum;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Fusion Artifacts (NABF) changed by B. K. Shreyamsha Kumar .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:p
   for jj=1:q
      if(gF(ii,jj)>gA(ii,jj) & gF(ii,jj)>gB(ii,jj))
         na(ii,jj)=1;
      else
         na(ii,jj)=0;
      end
   end
end
NABF=sum(sum(na.*((1-QAF).*wtA+(1-QBF).*wtB)))/wt_sum;

QABF+LABF+NABF1;
QABF+LABF+NABF;

author: Linfeng-Tang

repo: Image-Fusion


text:
function hole_point = p04_hole_point ( hole_index, m )

%*****************************************************************************80
%
%% P04_HOLE_POINT returns a point inside a given hole in problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer HOLE_INDEX, the index of the hole.
%
%    Input, integer M, the spatial dimension.
%
%    Output, real HOLE_POINT(M), a point in the hole
%
  hole_point(1:2) = [ 0.0, 0.0 ];

  return
end


text:
function headmodel = ft_headmodel_openmeeg(bnd, varargin)

% FT_HEADMODEL_OPENMEEG creates a volume conduction model of the
% head using the boundary element method (BEM). This function takes
% as input the triangulated surfaces that describe the boundaries and
% returns as output a volume conduction model which can be used to
% compute leadfields.
%
% This function implements
%   Gramfort et al. OpenMEEG: opensource software for quasistatic
%   bioelectromagnetics. Biomedical engineering online (2010) vol. 9 (1) pp. 45
%   http://www.biomedical-engineering-online.com/content/9/1/45
%   doi:10.1186/1475-925X-9-45
% and
%   Kybic et al. Generalized head models for MEG/EEG: boundary element method
%   beyond nested volumes. Phys. Med. Biol. (2006) vol. 51 pp. 1333-1346
%   doi:10.1088/0031-9155/51/5/021
%
% This link with FieldTrip is derived from the OpenMEEG project
% with contributions from Daniel Wong and Sarang Dalal, and uses external
% command-line executables. See http://openmeeg.github.io/
%
% Use as
%   headmodel = ft_headmodel_openmeeg(bnd, ...)
%
% Optional input arguments should be specified in key-value pairs and can
% include
%   conductivity     = vector, conductivity of each compartment
%   tissue           = tissue labels for each compartment
%
% See also FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD, FT_SYSMAT_OPENMEEG,
%          FT_SENSINTERP_OPENMEEG

%$Id$

ft_hastoolbox('openmeeg', 1);  % add to path (if not yet on path)
openmeeg_license;              % show the license (only once)
prefix = om_checkombin;        % check the installation of the binaries
if(~ispc) % if Linux/Mac, set number of threads
    omp_num_threads = feature('numCores');
    prefix = ['export OMP_NUM_THREADS=' num2str(omp_num_threads) ' && ' prefix];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the first part is largely shared with the dipoli and bemcp implementation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get the optional arguments
conductivity    = ft_getopt(varargin, 'conductivity');
tissue          = ft_getopt(varargin, 'tissue');

% copy the boundaries from the mesh into the volume conduction model
if isfield(bnd, 'bnd')
    bnd = bnd.bnd;
end

% rename pnt into pos
bnd = fixpos(bnd);

% start with an empty volume conductor
headmodel = [];

% determine the number of compartments
numboundaries = length(bnd);

% determine the desired nesting of the compartments
order = surface_nesting(bnd, 'outsidefirst');

% rearrange boundaries and conductivities
if numel(bnd)>1
    fprintf('reordering the boundaries to: ');
    fprintf('%d ', order);
    fprintf('\n');
    % update the order of the compartments
    bnd = bnd(order);
end

if isempty(conductivity)
    ft_warning('No conductivity is declared, using default values\n')
    if numboundaries == 1
        conductivity = 1;
    elseif numboundaries == 3
        % skin/skull/brain
        conductivity = [0.33 0.0042 0.33]
    elseif numboundaries == 4
        conductivity = [0.33 0.0042 1 0.33]
    else
        ft_error(['Conductivity values are required for ' num2str(numboundaries) ' shells'])
    end
    headmodel.cond = conductivity;
else
    if numel(conductivity)~=numboundaries
        ft_error('a conductivity value should be specified for each compartment');
    end
    % update the order of the compartments
    headmodel.cond = conductivity(order);
end

% assign default tissue labels if none provided
if(isempty(tissue))
    switch(numboundaries)
        case 3
            tissue = {'scalp','skull','brain'};
        case 4
            tissue = {'scalp','skull','CSF','brain'};
        otherwise
            tissue = strcat({'domain'},num2str((1:numboundaries)'));
    end
else
    tissue = tissue(order);
end
headmodel.tissue = tissue;

headmodel.skin_surface = 1;
headmodel.source = numboundaries;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this uses an implementation that was contributed by INRIA Odyssee Team
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

workdir = fullfile(tempdir,['ft_om_' datestr(now,'ddmmyyHHMMSSFFF')]);
mkdir(workdir);

try
    % Write the triangulations to file, named after tissue type.
    % OpenMEEG v2.3 and up internally adjusts the convention for surface
    % normals, but OpenMEEG v2.2 expects surface normals to point inwards;
    % this checks and corrects if needed
    bndfile = fullfile(workdir,strcat(tissue,'.tri'));
    for ii=1:length(bnd)
        ok = checknormals(bnd(ii));
        if ~ok
            % Flip faces for openmeeg convention (inwards normals)
            bndtmp = bnd(ii);
            bndtmp.tri = fliplr(bnd(ii).tri);

            ok = checknormals(bndtmp);
            if(ok)
                fprintf('flipping normals for OpenMEEG''s convention\n')
                bnd(ii).tri = bndtmp.tri;
            else
                % if still not ok after flip, throw a warning
                ft_warning('neither orientation of surface normals passes check... leaving as is.')
            end
            clear bndtmp
        end

        om_save_tri(bndfile{ii}, bnd(ii).pos, bnd(ii).tri);
    end

    % retain surfaces in headmodel structure (after possible normal flip)
    headmodel.bnd = bnd;

    condfile  = fullfile(workdir,'om.cond');
    geomfile  = fullfile(workdir,'om.geom');
    hmfile    = fullfile(workdir,'hm.bin');
    hminvfile = fullfile(workdir,'hminv.bin');

    % write conductivity and mesh files
    bndlabel = {};
    for i=1:length(bnd)
        [dum,bndlabel{i}] = fileparts(bndfile{i});
    end

    om_write_geom(geomfile,bndfile,bndlabel);
    om_write_cond(condfile,headmodel.cond,bndlabel);

    om_status = system([prefix 'om_assemble -HM ' geomfile ' ' condfile ' ' hmfile]);
    if(om_status ~= 0) % status = 0 if successful
        ft_error('Aborting OpenMEEG pipeline due to above error.');
    end

    headmodel.mat = inv(om_load_sym(hmfile,'binary'));

    rmdir(workdir,'s'); % remove workdir with intermediate files
catch
    disp(lasterr);
    rmdir(workdir,'s'); % remove workdir with intermediate files
    ft_error('an error occurred while running OpenMEEG');
end

% remember the type of volume conduction model
headmodel.type = 'openmeeg';

function ok = checknormals(bnd)
points = bnd.pos;
faces = bnd.tri;

% FIXME: this method is rigorous only for star shaped surfaces
ok = 0;
% translate to the center
org = mean(points,1);
points(:,1) = points(:,1) - org(1);
points(:,2) = points(:,2) - org(2);
points(:,3) = points(:,3) - org(3);

w = sum(solid_angle(points, faces));

if w<0 && (abs(w)-4*pi)<1000*eps
    ok = 0;
    disp('Your surface normals are outwards oriented...')
elseif w>0 && (abs(w)-4*pi)<1000*eps
    ok = 1;
    %   ft_warning('your normals are inwards oriented')
else
    ft_error('your surface probably is irregular\n')
    ok = 0;
end

author: spm


text:
function []=ps_scn_filt()
%PS_SCN_FILT estimate spatially correlated noise in unwrapped phase
%
%   Andy Hooper, June 2006
%
%   =======================================================================
%   11/2006 AH: Error corrected that was leaving master in temporal smoothing
%   04/2007 AH: Added 64-bit machine compatibility
%   05/2007 AH: Spatially correlated look angle error added
%   02/2010 AH: Replace unwrap_ifg_index with drop_ifg_index
%   02/2010 AH: Bug fixed that could cause a slave scn to be set to zero
%   11/2010 AH: If ramps estimated in step 7, subtract before scn estimation
%   02/2011 DB: Decreased time required for spatial filtering by factor 10
%   =======================================================================
logit;
fprintf('Estimating other spatially-correlated noise...\n')

pix_size=getparm('unwrap_grid_size',1);
time_win=getparm('scn_time_win',1);
deramp_ifg=getparm('scn_deramp_ifg',1);
scn_wavelength=getparm('scn_wavelength',1);
drop_ifg_index=getparm('drop_ifg_index',1);
small_baseline_flag=getparm('small_baseline_flag',1);

load psver
psname=['ps',num2str(psver)];
phuwname=['phuw',num2str(psver)];
sclaname=['scla',num2str(psver)];
apsname=['aps',num2str(psver)];
scnname=['scn',num2str(psver)]; % spatially-correlated noise

ps=load(psname);
uw=load(phuwname);
%aps=load(apsname);

if strcmpi(small_baseline_flag,'y')
    unwrap_ifg_index=[1:ps.n_image];
else
    unwrap_ifg_index=setdiff([1:ps.n_ifg],drop_ifg_index);
end

day=ps.day(unwrap_ifg_index);
master_ix=sum(ps.master_day>ps.day)+1;
n_ifg=length(unwrap_ifg_index);
n_ps=ps.n_ps;

ph_all=single(uw.ph_uw(:,unwrap_ifg_index));
if exist([sclaname,'.mat'],'file')
    scla=load(sclaname);
    ph_all=ph_all-single(scla.ph_scla(:,unwrap_ifg_index));
    ph_all=ph_all-repmat(single(scla.C_ps_uw),1,length(unwrap_ifg_index));
    if ~isempty(scla.ph_ramp)
        ph_all=ph_all-single(scla.ph_ramp(:,unwrap_ifg_index));
    end
end
ph_all(isnan(ph_all))=0;

disp(sprintf('   Number of points per ifg: %d',n_ps))

nodename=['scnfilt.1.node'];
fid=fopen(nodename,'w');
fprintf(fid,'%d 2 0 0\n',n_ps);

for i=1:n_ps
    fprintf(fid,'%d %f %f\n',i,ps.xy(i,2),ps.xy(i,3));
end

fclose(fid);

system('triangle -e scnfilt.1.node > triangle_scn.log');

fid=fopen('scnfilt.2.edge','r');
header=str2num(fgetl(fid));
N=header(1);
edges_nz=zeros(N,4);
for i=1:N
    edges_nz(i,:)=str2num(fgetl(fid));
end
fclose(fid);

%%% deramp end ifgs (unlike aps, orbit errors not so random and end
%%% orbit errors can pass through the low-pass filter
if strcmpi(deramp_ifg,'all')
    deramp_ifg=1:ps.n_ifg;
end
deramp_ifg=intersect(deramp_ifg,unwrap_ifg_index);
deramp_ix=zeros(size(deramp_ifg));
ph_ramp=zeros(n_ps,length(deramp_ifg));

if ~isempty(deramp_ifg)
    fprintf('   deramping selected ifgs...\n')
    G=double([ones(n_ps,1),ps.xy(:,2),ps.xy(:,3)]);
    %G=double([ones(n_ps,1),ps.xy(:,2)]); % range only

    for i=1:length(deramp_ifg)
        i3=find(unwrap_ifg_index==deramp_ifg(i))
        deramp_ix(i)=i3;
        d=(ph_all(:,i3));
        m=G\double(d(:));
        ph_this_ramp=G*m;
        ph_all(:,i3)=ph_all(:,i3)-ph_this_ramp; % subtract ramp
        ph_ramp(:,i)=ph_this_ramp;
    end
    save(scnname,'ph_ramp')
end

%%% smooth in time using gaussian moving window
isnanix=isnan(uw.ph_uw);
uw.ph_uw(isnanix)=0;
dph=ph_all(edges_nz(:,3),:)-ph_all(edges_nz(:,2),:);
dph_lpt=zeros(size(dph));
n_edges=size(dph,1);

fprintf('   low-pass filtering pixel-pairs in time...\n')

for i1=1:n_ifg
    time_diff_sq=(day(i1)-day)'.^2;
    weight_factor=exp(-time_diff_sq/2/time_win^2);
    weight_factor(master_ix)=0; % leave out master
    weight_factor=weight_factor/sum(weight_factor);
    dph_lpt(:,i1)=sum(dph.*repmat(weight_factor,n_edges,1),2);
end

dph_hpt=dph-dph_lpt;  % leaves master APS - slave APS - slave noise (+ residue master noise)

ph_hpt=zeros(n_ps-1,n_ifg);
ref_ix=1;

A=sparse([[1:n_edges]';[1:n_edges]'],[edges_nz(:,2);edges_nz(:,3)],[-ones(n_edges,1);ones(n_edges,1)]);
A=double(A(:,[1:ref_ix-1,ref_ix+1:n_ps]));

fprintf('   solving for high-frequency (in time) pixel phase...\n')

for i=1:n_ifg
    ph_hpt(:,i)=A\double(dph_hpt(:,i));
end

ph_hpt=[ph_hpt(1:ref_ix-1,:);zeros(1,n_ifg);ph_hpt(ref_ix:end,:)]; % add back ref point

ph_hpt(:,deramp_ix)=ph_hpt(:,deramp_ix)+ph_ramp;

ph_hpt=single(ph_hpt);

sigma_sq_times_2=2*scn_wavelength.^2;
ph_scn=nan(n_ps,n_ifg);
patch_dist=scn_wavelength*4;
patch_dist_sq=patch_dist*patch_dist;
ix_range=ceil(n_ps/(max(ps.xy(:,3))-min(ps.xy(:,3)))*patch_dist*0.2);
ix1=1;
ix2=ix_range;
ps.xy(:,1)=[1:n_ps]';

fprintf('   low-pass filtering in space...\n')

for i=1:n_ps
    x_min=ps.xy(i,2)-patch_dist;
    x_max=ps.xy(i,2)+patch_dist;
    y_min=ps.xy(i,3)-patch_dist;
    y_max=ps.xy(i,3)+patch_dist;

    ix1=ix1+ix_range;
    ix1(ix1>n_ps)=n_ps;
    while ix1>1 & ps.xy(ix1-1,3)>=y_min
        ix1=ix1-ix_range;
    end

    ix2=ix2-ix_range;
    ix2(ix2<1)=1;
        while ix2<n_ps & ps.xy(ix2+1,3)<=y_max
        ix2=ix2+ix_range;
    end

    ix1(ix1<1)=1;
    ix2(ix2>n_ps)=n_ps;

    xy_near=ps.xy(ix1:ix2,:);
    xy_near=xy_near(xy_near(:,2)>=x_min & xy_near(:,2)<=x_max & xy_near(:,3)>=y_min & xy_near(:,3)<=y_max,:);
    dist_sq=(xy_near(:,2)-ps.xy(i,2)).^2+(xy_near(:,3)-ps.xy(i,3)).^2;
    in_range_ix=dist_sq<patch_dist_sq; % exclude those out of range
    xy_near=xy_near(in_range_ix);
    dist_sq=dist_sq(in_range_ix);
    weight_factor=exp(-dist_sq/sigma_sq_times_2);
    weight_factor=weight_factor/sum(weight_factor); % normalize

    %ph_scn(i,:)=sum(ph_hpt(xy_near(:,1),:).*repmat(weight_factor,1,n_ifg),1); 		% slow method
    ph_scn(i,:)=weight_factor'*ph_hpt(xy_near(:,1),:);  				% speed increased by a factor 10

    if i/1000==floor(i/1000)
        disp([num2str(i),' PS processed'])
    end
end

ph_scn=ph_scn-repmat(ph_scn(1,:),n_ps,1); % re-ref to 1st PS
ph_scn_slave=zeros(size(uw.ph_uw));
ph_scn_slave(:,unwrap_ifg_index)=ph_scn;
ph_scn_slave(:,master_ix)=0;

save(scnname,'ph_scn_slave','ph_hpt','ph_ramp')
logit(1);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function test_example_sourcemodel_aligned2mni_atlas

% MEM 4gb
% WALLTIME 00:10:00

%
%% Create brain atlas based MNI-aligned grids in individual head-space
%
% When combining the source-level data of multiple subjects this data is typically first interpolated (e.g., using **[ft_sourceinterpolate](https://github.com/fieldtrip/fieldtrip/blob/release/ft_sourceinterpolate.m)**) and then spatially normalized to a template brain (e.g., using **[ft_volumenormalise](https://github.com/fieldtrip/fieldtrip/blob/release/ft_volumenormalise.m)**). However it is also possible to define the source-reconstruction grid for each individual subject in such a way that all these grids are already aligned in MNI-space. The combination or statistic of source-level data across subjects can then directly be computed within the source-structure without the need to interpolate and normalize each volume. In addition, the position of the grid points can be chosen in a way that only locations corresponding to particular brain areas (parcels) are included.
%
% The procedure is as follows. First, a template grid is computed on the basis of the standard head model located in the FieldTrip template directory. Subsequently, a brain atlas is loaded using **[ft_read_atlas](https://github.com/fieldtrip/fieldtrip/blob/release/fileio/ft_read_atlas.m)**. On the basis on the brain atlas and the template grid it is possible to create a binary mask of all locations in the template grid that match atlas locations using **[ft_volumelookup](https://github.com/fieldtrip/fieldtrip/blob/release/ft_volumelookup.m)**. Finally the mask is used to determine which location shall be defined as 'inside' the head model.
%
[ftver, ftpath] = ft_version;
load(fullfile(ftpath, 'template/headmodel/standard_singleshell.mat'));

cfg = [];
cfg.xgrid  = -20:1:20;
cfg.ygrid  = -20:1:20;
cfg.zgrid  = -20:1:20;
cfg.unit   = 'cm';
cfg.tight  = 'yes';
cfg.inwardshift = -1.5;
cfg.headmodel   = vol;
template_grid   = ft_prepare_sourcemodel(cfg);
template_grid.coordsys = 'mni';

template_grid = ft_convert_units(template_grid,'cm');

figure;
hold on
ft_plot_mesh(template_grid.pos(template_grid.inside,:));
ft_plot_headmodel(vol,  'facecolor', 'cortex', 'edgecolor', 'none');
ft_plot_axes(vol);
alpha 0.5
camlight

% Read the atlas, convert to units of cm and create the binary mask.
%
atlas = ft_read_atlas(fullfile(ftpath, 'template/atlas/aal/ROI_MNI_V4.nii'));

atlas = ft_convert_units(atlas,'cm');

cfg = [];
cfg.atlas      = atlas;
cfg.roi        = atlas.tissuelabel;  % here you can also specify a single label, i.e. single ROI
mask           = ft_volumelookup(cfg, template_grid);

% Now we determine all indices of the binary mask to be considered as inside the head model. And plot the result. Note the missing dipole locations for example in the vicinity of the ventricles.
%
template_grid.inside = false(template_grid.dim);
template_grid.inside(mask==1) = true;

figure;
ft_plot_mesh(template_grid.pos(template_grid.inside,:));

%
% Load the subject-specific MRI from [here](ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/salzburg/mri.mat) and inverse-warp the subject specific grid to the template grid.
%
mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/salzburg/mri.mat'));

cfg                = [];
cfg.warpmni   = 'yes';
cfg.template  = template_grid;
cfg.nonlinear = 'yes'; % use non-linear normalization
cfg.mri            = mri;
sourcemodel        = ft_prepare_sourcemodel(cfg);

% Finally, you can load the subject-specific headmodel from [here](ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/salzburg/hdm.mat) and check the result with the following code.
%
close all

hdm = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/salzburg/hdm.mat'));

hdm         = ft_convert_units(hdm, 'm');
sourcemodel = ft_convert_units(sourcemodel, 'm');

figure
hold on
ft_plot_headmodel(hdm,  'facecolor', 'cortex', 'edgecolor', 'none');
ft_plot_axes(hdm);
alpha 0.4  % make the surface transparent
ft_plot_mesh(sourcemodel.pos(sourcemodel.inside,:));
view ([0 -90 0])

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef refractiveIndexTensor < tensor

  methods

    function rI = refractiveIndexTensor(varargin)
      rI = rI@tensor(varargin{:},'rank',2);

      % TODO: set up the unit correctly !!!

    end

  end

  methods (Static = true)
    function rI = calcite
      cs = crystalSymmetry('-3m1',[5,5,17],'mineral','Calcite','X||a');
      rI = refractiveIndexTensor(diag([1.66 1.66 1.486]),cs);
    end

    function rI = olivin
      cs = loadCIF('olivin.cif');
      rI = refractiveIndexTensor(diag([1.640 1.660 1.680]),cs);
    end

    function test
      rI = refractiveIndexTensor.calcite;

      vprop = plotS2Grid;

      figure(1)
      n = rI.birefringence(vprop);

      plot3d(vprop,n)

      figure(2)
      thickness = 10000;
      rgb = spectralTransmission(rI,vprop,thickness);
      plot3d(vprop,rgb./100)

    end

    function test2
      mtexdata fo

      rI = refractiveIndexTensor.olivin;

      oM = spectralTransmissionColorKey(rI,1000);
      oM.polarizer = vector3d.Y;

      plot(ebsd('fo'),oM.orientation2color(ebsd('fo').orientations)./100)

    end

  end

end

author: mtex-toolbox

repo: mtex

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


text:
function [ n_data, alpha, beta, x, fx ] = extreme_values_cdf_values ( n_data )

%*****************************************************************************80
%
%% EXTREME_VALUES_CDF_VALUES returns some values of the Extreme Values CDF.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Needs["Statistics`ContinuousDistributions`"]
%      dist = ExtremeValuesDistribution [ alpha, beta ]
%      CDF [ dist, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ALPHA, the first parameter of the distribution.
%
%    Output, real BETA, the second parameter of the distribution.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 12;

  alpha_vec = [ ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.4000000000000000E+01, ...
     0.5000000000000000E+01 ];

  beta_vec = [ ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.4000000000000000E+01, ...
     0.5000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01 ];

  fx_vec = [ ...
     0.3678794411714423E+00, ...
     0.8734230184931166E+00, ...
     0.9818510730616665E+00, ...
     0.9975243173927525E+00, ...
     0.5452392118926051E+00, ...
     0.4884435800065159E+00, ...
     0.4589560693076638E+00, ...
     0.4409910259429826E+00, ...
     0.5452392118926051E+00, ...
     0.3678794411714423E+00, ...
     0.1922956455479649E+00, ...
     0.6598803584531254E-01 ];

  x_vec = [ ...
     0.1000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.4000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.3000000000000000E+01 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

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

  return
end


text:
function pols = kjacopols ( x, a, b, n )

%*****************************************************************************80
%
%% KJACOPOLS evaluates Jacobi polynomials.
%
%  Discussion:
%
%    This routine evaluates Jacobi polynomials.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    26 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Input, real A, B, the parameter values.
%
%    Input, integer N, the highest degree to be evaluated.
%
%    Output, real POLS(N+1), the polynomial values.
%
  pols = zeros(n+1,1);

  pkp1 = 1.0;
  pols(1) = pkp1;

  if ( n == 0 )
    return
  end

  pk = pkp1;
  pkp1 = ( a / 2.0 - b / 2.0 ) ...
    + ( 1.0 + a / 2.0 + b / 2.0 ) * x;
  pols(2) = pkp1;

  if ( n == 1 )
    return
  end

  for k = 2 : n

    pkm1 = pk;
    pk = pkp1;

    alpha = ( 2.0 * k + a + b - 1.0 ) ...
      * ( a * a - b * b + ( 2.0 * k + a + b - 2.0 ) ...
      * ( 2.0 * k + a + b ) * x );

    beta = 2.0 * ( k + a - 1.0 ) * ( k + b - 1.0 ) ...
      * ( 2.0 * k + a + b );

    pkp1 = ( alpha * pk - beta * pkm1 ) ...
      / ( 2.0 * k * ( k + a + b ) ...
      * ( 2.0 * k + a + b - 2.0 ) );

    pols(k+1) = pkp1;

  end

  return
end


text:
%% IPOPT PARFOR testing
clc
clear all

%Objective
fun = @(x) 20 + x(1)^2 + x(2)^2 - 10*(cos(2*pi*x(1)) + cos(2*pi*x(2)));
%Constraints
lb = [5*pi;-20*pi];
ub = [20*pi;-4*pi];
%Setup Options
opts = optiset('solver','ipopt');
Opt = opti('fun',fun,'bounds',lb,ub,'opts',opts);

%Generate a series of starting points
n = 1000;
X0 = [linspace(lb(1),lb(end),n); linspace(ub(1),ub(end),n)];

%Preallocate solution vector
sols = zeros(2,n);
solp = zeros(2,n);

%% run serial (8.17s)
tic
for i = 1:n
    sols(:,i) = solve(Opt,X0(:,i));
end
toc

%%
matlabpool(4);

%% run parfor (2.6s)
tic
parfor i = 1:n
    solp(:,i) = solve(Opt,X0(:,i));
end
toc

%% check results
err = norm(sols-solp)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% PlotNormScatt(x,Scatt)
% Color plot of a normalized scattering representation
% x gives the horizontal axis
% nScatt is the normalized scattering tranform which is piecewise constant
% The color on each interval depends upon the path length:
% 0 is yellow, 1 is red, 2 is green 3 is blue, 4 is magenta

function plot_diracnorm_scat(x,nScatt)
	color(1) = 'y';
	color(2) = 'r';
	color(3) = 'g';
	color(4) = 'b';
	color(5) = 'm';
	color(6) = 'k';
	color(7:99) = 'k';

	N=size(x,2);

	n=1;
	while(n<N)
		X(1) = x(n);
		Y(1) = nScatt(n);
		m = (round(10^(7) * Y(1) - 10 * floor(Y(1) * 10^(6))))+1;
		y = Y(1);
		while ((y == Y(1)) && (n < N))
			n = n+1;
			y = nScatt(n);
		end
		X(2) = x(n-1);
		Y(2) = Y(1);
		plot(X,Y,color(m));
		hold on;
		X(1) = X(2);
		Y(1) = y;
		plot(X,Y,color(m));
	end
	hold off;
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [xir,xjc] = sdpasplit(xcol,maxnnz)
% ------------------------------------------------------------
% Allocate enough space
% ------------------------------------------------------------
 n = length(xcol);
 if nargin < 2
   maxnnz = n;
 end
 xir = zeros(maxnnz,1);
 xjc = zeros(maxnnz+1,1);
% ------------------------------------------------------------
% Find all nonempty blocks, and store their block start in xjc.
% ------------------------------------------------------------
 knz = 0;
 k = -1;
 for j = 1:n
   if xcol(j) > k
     knz = knz + 1;
     k = xcol(j);
     xir(knz) = k;
     xjc(knz) = j;
   end
 end
% ------------------------------------------------------------
% Close last block in xjc and shrink xjc, xir to correct size
% ------------------------------------------------------------
 xjc(knz+1) = n+1;
 xjc = xjc(1:knz+1);
 xir = xir(1:knz);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
figure(1)
plot(t, Y{3}')

author: PatternRecognition

repo: OpenBMI


text:
% Discrete Cepstral Envelope (DCE)
%
% Input
%  af            : Array of structures with matrices containing sinusoidal
%                  parameters (as the output of sin_analysis.m).
%                  Each matrix is made of:
%                   The 1st line for the frequency of each sinusoid [Hz]
%                   The 2nd line for their amplitude (linear scale)
%                  The DC has to be included (in the first column)
%  fs            : [Hz] Signal's sampling frequency
%  order         : Cepstral order
%  [extrap_dcny] : If true (default), alleviate stability problems of the DCE
%                     by replacing the DCE value and extrapolating sinusoidal
%                     components up to Nyquist.
%                  If false, use the sinusoidal components as they are.
%  [scale]       : empty  : Frequency linear scale (default)
%                 'mel'  : see frq2mel (for MFCC computation)
%                 'bark' : see frq2bark
%                 'erb'  : see frq2erb
%  [Bw]          : [Hz] Standard-deviation of the Gaussian used as weighting
%                       function (for emphasizing the importance of the low
%                       frequencies in the solution).
%                       Bw As to be big enough to have the weights still
%                       significant up to Nyquist.
%  [lr]          : Regularization parameter (as in [2]) (def. 0)
%  [dftlen]      : DFT's length, if the 4th output argument is requested.
%
% Output
%  cc             : Cepstral coefficients
%  E              : The amplitude cepstral envelope
%
% References
%  [1] T. Galas and X. Rodet, "Generalized discrete cepstral analysis for
%      deconvolution of source-filter system with discrete spectra,"
%      in IEEE Applications of Signal Processing to Audio and Acoustics (ASSP)
%      Workshop, pp. 71-72, 1991.
%  [2] M. Campedel-Oudot, O. Cappe and E. Moulines, "Estimation of the Spectral
%      Envelope of Voiced Sounds Using a Penalized Likelihood Approach", IEEE
%      Transactions on Speech and Audio Processing, vol. 9, no. 5, pp. 469-481,
%      2001.
%
% Copyright (c) 2011 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the COVAREP project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [cc E] = env_dce_sfa(af, fs, order, extrap_dcny, scale, Bw, lr, dftlen)

    % Input parameters
    if nargin<4; extrap_dcny=true; end
    if nargin<5; scale = []; end
    if nargin<6 || isempty(Bw); Bw = 2000*(fs/16000); end % To make it as
                                           % similar as possible to the DCE-MFA
    if nargin<7 || isempty(lr); lr = 0.035; end % [2]
    if nargin<8; dftlen=4096; end

    if ~isempty(scale)
        Bw = 100*fs;
        lr = 0.035; % To make it as similar as possible to [2]
        eval(['fnscale=@frq2' scale ';']);
        af.f = 0.5*fs*fnscale(af(n).f)/fnscale(fs/2);
    end

    % If asked, extrapolate components at DC and up to Nyquist
    if extrap_dcny; af = env_extrap_sins_dcny(af, fs); end

    af.f = af.sins(1,:);
    af.a = log(af.sins(2,:));

    % Prepare the weighting functions
    fk = af.f(:);
    Nk = length(fk);
    wk = exp(-fk.^2 / (2 * Bw * Bw))';
    wk = wk./sum(wk);
    Wk = diag(wk);

    % Compute the DCE envelope
    fk = af.f(:);
    ak = af.a(:);
    Nk = length(fk);
    Bk = [ones(Nk,1) cos(2*pi*fk/fs*(1:order))];

    BWB = (Bk'*Wk*Bk);
    rt =  (Bk'*Wk*(ak));

    if lr>0
        cc = (BWB+lr*diag(ones(size(BWB,1),1)))\rt;
    else
        cc = BWB\rt;
    end

    % If asked, compute the envelope
    if nargout>1
        if isempty(scale)
            E = fft(cc, dftlen);
            E = exp(E(1:end/2+1));
        else
            if strcmp(scale,'bark')
                E = barkcc2spec(cc, fs, dftlen);
                E = E(1:end/2+1);
            elseif strcmp(scale,'mel')
                E = exp(fft(cc, dftlen));
                E = E(1:dftlen/2+1);
                E = cc2hspec(cc, fs, dftlen);
                E = fwcep2hspec(cc, fs, dftlen);
            end
        end
    end

    % Plot final solution
    if 0
        hold off;

        global S Erefg selfr;
        if ~isempty(S);
            Sbins = fs*(0:dftlen/2)/dftlen;
            plot(Sbins, mag2db(abs(S)), ':k');
            hold on;
        else; hold off; end
        if ~isempty(Erefg);
            Erefgbins = 0.5*fs*(0:length(Erefg)-1)/length(Erefg);
            plot(Erefgbins, mag2db(abs(Erefg)), 'k');
            hold on;
        end

        plot(af.f, mag2db(exp(af.a)),'xr');
        hold on;
        F = fs*(0:dftlen/2)/dftlen;
%          plot(F, mag2db(abs(L)), 'g');
        hold on;
        plot(F, mag2db(abs(E)), 'b', 'LineWidth', 2);
        keyboard
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%% MEEG time generalization multiple comparison correction
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm where a participant saw
% images of six object categories.
%
% The code presented here can be adapted for other MEEG analyses, but
% there please note:
% * the current examples do not perform baseline corrections or signal
%   normalizations, which may reduce discriminatory power.
%
% Note: running this code requires FieldTrip.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% get timelock data in CoSMoMVPA format

% set configuration
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'meg_obj6');

% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);

% reset citation list
cosmo_check_external('-tic');

% load preprocessed data
data_fn=fullfile(data_path,'meg_obj6_s00.mat');
data_tl=load(data_fn);

% convert to cosmomvpa struct and show the dataset
ds=cosmo_meeg_dataset(data_tl);
cosmo_disp(ds);

% set the targets (trial condition)
ds.sa.targets=ds.sa.trialinfo(:,1); % 6 categories

% set the chunks (independent measurements)
% all trials are here considered to be independent
nsamples=size(ds.samples,1);
ds.sa.chunks=(1:nsamples)';

% in addition give a label to each trial
index2label={'body','car','face','flower','insect','scene'};
ds.sa.labels=cellfun(@(x)index2label(x),num2cell(ds.sa.targets));

% just to check everything is ok
cosmo_check_dataset(ds);

%% Select subset of sensors and time points

% Select posterior gradiometers
sensor_posterior_planar={'MEG1632', 'MEG1642', 'MEG1732', 'MEG1842', ...
                        'MEG1912', 'MEG1922', 'MEG1942', 'MEG2232', ...
                        'MEG2312', 'MEG2322', 'MEG2342', 'MEG2432', ...
                        'MEG2442', 'MEG2512', 'MEG2532',...
                        'MEG1633', 'MEG1643', 'MEG1733', 'MEG1843', ...
                        'MEG1913', 'MEG1923', 'MEG1943', 'MEG2233', ...
                        'MEG2313', 'MEG2323', 'MEG2343', 'MEG2433', ...
                        'MEG2443', 'MEG2513', 'MEG2533'};

msk=cosmo_dim_match(ds,'chan',sensor_posterior_planar,...
                        'time',@(t)t>=0 & t<=.3);

ds_sel=cosmo_slice(ds,msk,2);
ds_sel=cosmo_dim_prune(ds_sel);

%%
% subsample time dimension to speed up the analysis
subsample_time_factor=3; % take every 3rd time point

% take every subsample_time_factor-th
% hint: use ds.fa.time to find the desired features, then use cosmo_slice
% and cosmo_dim_prune.
% >@@>
msk = mod(ds_sel.fa.time,subsample_time_factor)==1;
ds_sel=cosmo_slice(ds_sel,msk,2);
ds_sel=cosmo_dim_prune(ds_sel);
% <@@<

% to illustrate group analysis, we use data from a single participant
% and divide it in ten parts. Each part represents a pseudo-participant.
n_pseudo_participants = 10;
ds_sel.sa.subject_id=cosmo_chunkize(ds_sel,n_pseudo_participants);
ds_cell=cosmo_split(ds_sel, 'subject_id');

%%
% apply the cosmo_dim_generalization_measure to the data from each
% pseudo-participant
group_cell=cell(n_pseudo_participants,1);

for k=1:n_pseudo_participants
    ds_subj=ds_cell{k};
    ds_subj.sa=rmfield(ds_subj.sa,'subject_id');

    ds_subj=cosmo_balance_dataset(ds_subj);
    ds_subj.sa.chunks=cosmo_chunkize(ds_subj,2);
    ds_subj_tr=cosmo_dim_transpose(ds_subj,'time',1);

    % use a custom measure that computes a one-way ANOVA F-value and
    % then converts this to a z-score
    measure=@(d,opt)cosmo_stat(d,'F','z');

    ds_time_gen=cosmo_dim_generalization_measure(ds_subj_tr,...
                                'measure',@cosmo_correlation_measure,...
                                'dimension','time');

    group_cell{k} = ds_time_gen;
end

%%
% show an element of group_cell. What is the size of .samples?
cosmo_disp(group_cell{1});

%%
% To do group analysis, the above format will not work.
% We want a dataset ds_group with size n_pseudo_participants x NF
% where NF is the number of features.

% allocate a cell group_cell_tr with the same size of group_cell
% >@@>
group_cell_tr=cell(size(group_cell));
% <@@<

for k=1:numel(group_cell)
    % take data from the k-th participant and store
    % in a varibale ds_time_gen
    ds_time_gen=group_cell{k};

    % change 'train_time' and 'test_time' from being sample dimensions
    % to become feature dimensions.
    % Hint: use cosmo_dim_transpose.
    % >@@>
    ds_time_gen_tr=cosmo_dim_transpose(ds_time_gen,...
                                        {'train_time','test_time'},2);
    % <@@<

    % set chunks and targets for a one-sample t-test against zero,
    % so that across participants: all targets have the same value, and
    % all chunks have different values.
    % >@@>
    ds_time_gen_tr.sa.chunks=k;
    ds_time_gen_tr.sa.targets=1;
    % <@@<

    % store ds_time_gen_tr as the k-th element in group_cell_tr
    % >@@>
    group_cell_tr{k}=ds_time_gen_tr;
    % <@@<
end

% show an element of group_cell_tr. What is the size of .samples?
cosmo_disp(group_cell_tr{1});

%%
% stack the elements in group_cell_tr into a dataset ds_group
% >@@>
ds_group=cosmo_stack(group_cell_tr);
% <@@<

%%
% define a clustering neighborhood and store the result in a struct
% called nbrhood
% Hint: use cosmo_cluster_neighborhood
nbrhood=cosmo_cluster_neighborhood(ds_group);

%%

% run multiple comparison correction using cosmo_montecarlo_cluster_stat
% with 1000 iterations, for a t-test against h0_mean=0.
opt=struct();
opt.niter=1000;
opt.h0_mean=0;
ds_tfce=cosmo_montecarlo_cluster_stat(ds_group,nbrhood,opt);

%%
% >@@>
ds_group=cosmo_stack(group_cell_tr);
% <@@<

% extract the values from the ds_tfce, using cosmo_unflatten.
% store the array, and the dimension labels and values, into variables
% arr, dim_labels, and dim_values.
% Hint: use cosmo_unflatten.
% >@@>
[arr, dim_labels, dim_values]=cosmo_unflatten(ds_tfce);
% <@@<

% Reshape to arr to be 2-dimensional and store the result in arr_2d,
% then visualize the array
% >@@>
arr_2d=squeeze(arr);
% <@@<
clim=[-1,1]*max(abs(arr_2d(:)));
imagesc(arr_2d,clim);

% add axis labels
nticks=5;

ytick=round(linspace(1, numel(dim_values{1}), nticks));
ylabel(strrep(dim_labels{1},'_',' '));
set(gca,'Ytick',ytick,'YTickLabel',dim_values{1}(ytick));

xtick=round(linspace(1, numel(dim_values{2}), nticks));
xlabel(strrep(dim_labels{2},'_',' '));
set(gca,'Xtick',xtick,'XTickLabel',dim_values{2}(xtick));
colorbar();

% bonus: add markers indicating significance
% >@@>
z_min=1.96;
[i,j]=find(abs(arr_2d)>z_min);
hold on;
scatter(j,i,'o','k');
hold off;
% <@@<

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [y,zf]=ditherq(x,m,zi)
%DITHERQ  add dither and quantize [Y,ZF]=(X,M,ZI)
%  Inputs:
%      x   is the input signal
%	   m   specifies the mode:
%          'w'  white dither (default)
%          'h'  high-pass dither (filtered by 1 - z^-1)
%          'l'  low pass filter  (filtered by 1 + z^-1)
%          'n'  no dither

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: ditherq.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=size(x);
n=length(x);
if nargin<3 | ~length(zi)
    zi=rand(1);
end
    if nargin<2
        m='w';
    end
if any(m=='n')
    y=round(x);
elseif any(m=='h') | any(m=='l')
    v=rand(n+1,1);
    v(1)=zi;
    zf=v(end);
    if any(m=='h')
        y=round(x(:)+v(2:end)-v(1:end-1));
    else
        y=round(x(:)+v(2:end)+v(1:end-1)-1);
    end
else
    y=round(x(:)+rand(n,2)*[1;-1]);
    zf=rand(1);                         % output a random number anyway
end
if s(1)==1
    y=y.';
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
%% ------------------------solveRWF.m--------------------------------------

% Solver for Reweighted Wirtinger Flow as given in Algorithm 1
% of the Reweightde Wirtinger Flow (TAF) paper. Refer to the userguide for
% a detailed usage of the package.

% PAPER TITLE:
%              Phase Retrieval via Reweighted Wirtinger Flow
% ARXIV LINK:
%              https://doi.org/10.1364/AO.56.002418

% INPUTS:
%         A:   Function handle/numerical matrix for data matrix A. The rows
%              of this matrix are the measurement vectors that produce
%              amplitude measurements '\psi'.
%         At:  Function handle/numerical matrix for A transpose.
%         b0:  Observed data vector consisting of amplitude measurements
%              generated from b0 = |A*x|. We assign it to 'psi' to be
%              consistent with the notation in the paper.
%         x0:  The initial vector to be used by any solver.
%        opts: struct consists of the options for the algorithm. For
%              details,see header in solvePhaseRetrieval.m or the User
%              Guide.

% OUPTUT :
%         sol: n x 1 vector. It is the estimated signal.
%        outs: A struct consists of the convergence info. For details,
%              see header in solvePhaseRetrieval.m or the User Guide.

% Note:        When a function handle is used, the value of 'n' (the length
%              of the unknown signal) and 'At' (a function handle for the
%              adjoint of 'A') must be supplied. When 'A' is numeric, the
%              values of 'At' and 'n' are ignored and inferred from the
%              arguments

% DESCRIPTION:
%             The spectral method proposed by Candes et al. in their
%             seminal Wirtinger Flow paper suffers from fat tails that can
%             distort the true solution. A possible way to overcome this is
%             to truncate the outliers. In this paper, the authors propose
%             a weighted truncated iterative descent method that removes
%             those measurement vector a_m that are not correlated to the
%             intial guess. This truncation and weighting step is done at
%             each iteration of the gradient descent method.

% METHOD:
%         1) Our implementation uses FASTA, a fast gradient solver.
%
%         2) Set the objectve f = @(z) 1/2 * sum(weights .*
%             (abs(z) - b0).^2). The weights are computed according to
%             equation (14) in algorithm 1 of the proposing paper
%
%         3) The gradient is grad = @(z) weights .* (z - b0 .* sign(z)).
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START-----------------------------------

function [sol, outs] = solveRWF(A, At, b0, x0, opts)
%     addpath('solvers/linesearch');

    m = length(b0);

    innerOpts = struct;
    innerOpts.maxIters = opts.maxIters;
    innerOpts.maxTime = opts.maxTime;
    innerOpts.tol = opts.tol;
    innerOpts.verbose = opts.verbose;
    innerOpts.recordTimes = opts.recordTimes;
    innerOpts.recordResiduals = opts.recordResiduals;
    innerOpts.recordMeasurementErrors = opts.recordMeasurementErrors;
    innerOpts.recordReconErrors = opts.recordReconErrors;
    innerOpts.xt = opts.xt;

    innerOpts.updateObjectivePeriod = opts.reweightPeriod;
    innerOpts.searchMethod = opts.searchMethod;
    innerOpts.betaChoice = opts.betaChoice;

    [sol, outs] = gradientDescentSolver(A, At, x0, b0, @updateObjective, innerOpts);

    function [f, gradf] = updateObjective(~, Ax)
        weights = 1 ./ (abs(abs(Ax).^2 - b0.^2) + opts.eta*ones(m, 1));
        s = sum(weights);
        f = @(z) 0.5/s * sum(weights .* (abs(z).^2 - b0.^2).^2);
        gradf = @(z) (1.0/s)*weights .* (abs(z).^2 - b0.^2) .* z;
    end
end

author: tomgoldstein

repo: phasepack-matlab


text:
function json=get_coverage_json(obj, root_dir)
% Get JSON coverage representation
%
% json=get_coverage_json(obj, root_dir)
%
% Inputs:
%   obj                 MOcovMFile instance
%   root_dir            git root directory in which the file represented
%                       by obj resides
%
% Output:
%   json                JSON String representation of coverage, with
%                       elements:
%                       'name'          filename, relative to root_dir
%                       'source_digest' MD5 checksum of file contents
%                       'coverage'      Array with number of times each
%                                       line was executed; or null for
%                                       lines that cannot not be executed.
%
% Notes:
%   - this output can be used by the coveralls.io online coverage service
%     in combination with travis-ci

    name=mocov_get_relative_path(root_dir,obj.filename);
    source_digest=mocov_util_md5(obj.filename);
    coverage=get_coverage(obj);

    json=sprintf(['{ "name": "%s",\n',...
                    '"source_digest": "%s",\n'...
                    '"coverage": %s\n }\n'],...
                    name,source_digest,coverage);

function json_coverage=get_coverage(obj)
    executable=get_lines_executable(obj);
    executed_count=get_lines_executed_count(obj);

    n=numel(executable);

    json_parts=cell(1,2*n+1);
    json_parts{1}='[';
    for k=1:n
        line_executable=executable(k);

        if line_executable
            coverage_str=sprintf('%d',executed_count(k));
        else
            coverage_str='null';
        end

        json_parts{2*k}=coverage_str;
        if k<n
            suffix=',';
        else
            suffix=']';
        end

        json_parts{2*k+1}=suffix;
    end

    json_coverage=sprintf('%s',json_parts{:});

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function par = fmri_ldpar(varargin)
%
% Load specified par files for multiple runs.
%
% par = fmri_ldpar(ParFileList)
% par = fmri_ldpar(ParFile1, ParFile2, ...)
%
% ParFileList is a vertical cat of parfile names (ie, each run's'
% parfile name on a different row).
%
% ParFile format is as follows:
%   Column 1: stimulus presentation time (float or int)
%   Column 2: stimulus condition number
%   Column > 2: ignored
%
% par dimensionality: nPresentations x 2 x nRuns
%
% ras, 04/05: sorts rows by onsets, so they need not
% be specified in chronological order.
% ras, 02/06: replaced w/ mrVista tools, to fix linux \n issue.
par = [];

if(nargin == 0)
    error('USAGE: par = fmri_ldpar(ParFile)');
end

if( length(varargin) == 1)
    ParFile = varargin{1};
    nRuns = size(ParFile,1);
else
    nRuns = length(varargin);
    ParFile = '';
    for r = 1:nRuns,
        ParFile = strvcat(ParFile,varargin{r});
    end
end

% Go through each run and load par %
for r=1:nRuns,
    [onsets conds] = readParFile(deblank(ParFile(r,:)));
    par = [par; onsets(:) conds(:)];
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%--- help for statespace/forecast ---
%
% FORECAST Forecast states and observations of state-space models
%
%  Syntax:
%
%    [Y,YMSE,X,XMSE] = forecast(Mdl,numPeriods,Y0)
%    [Y,YMSE,X,XMSE] = forecast(Mdl,numPeriods,name,value,...)
%
%  Description:
%
%    Generate multiple-period forecasts for observation vector y(t) and state
%    vector x(t) for a general state-space model (SSM) of the form:
%
%    State equation:       x(t) = A(t) * x(t-1) + B(t) * u(t)
%    Observation equation: y(t) = C(t) * x(t)   + D(t) * e(t)
%
%    where u(t) and e(t) are uncorrelated, unit-variance white noise vector
%    processes. The length of x(t), y(t), u(t), and e(t) is m, n, k, and h,
%    respectively.
%
%    Additionally, forecast uncertainties are also generated.
%
%  Input Arguments:
%
%    Mdl - A state-space model, as created by the SSM constructor or
%      SSM/ESTIMATE method.
%
%    numPeriods - Positive, scalar, integer specifying the forecast horizon.
%
%    Y0 - Observed response data to be forecasted. For time-invariant models
%      in which the length of each observation vector (n) is the same, Y0 is
%      a T-by-n matrix. For time-varying models in which the length of the
%      observation vector changes, Y0 is a T-by-1 cell array in which each
%      element contains a time-varying n-element vector of observations, y(t),
%      associated with the corresponding period. The last observation is the
%      most recent.
%
%  Optional Input Name/Value Pairs:
%
%    'A' Cell vector of forecasted state transition matrices in which each
%        element is a matrix corresponding to a period in the forecast horizon
%        at time t = 1,2,...,numPeriods. If the length of the state vector
%        x(t) is constant, then each element is a square m-by-m matrix;
%        however, if the length of x(t) changes, then some elements are
%        non-square matrices. The length of the cell array must be at least
%        numPeriods, and any elements beyond the forecast horizon are ignored.
%        By default, the last coefficient of the input SSM model Mdl is used
%        for all future periods.
%
%    'B' Cell vector of forecasted state disturbance loading matrices in which
%        each element is a matrix corresponding to a period in the forecast
%        horizon at time t = 1,2,...,numPeriods. If the lengths of the state
%        vector x(t) and disturbance vector u(t) are constant, then each
%        element is an m-by-k matrix; however, if the length of x(t) or u(t)
%        changes, then the elements are matrices of various sizes. The length
%        of the cell array must be at least numPeriods, and any elements
%        beyond the forecast horizon are ignored. By default, the last
%        coefficient of the input SSM model Mdl is used for all future periods.
%
%    'C' Cell vector of forecasted measurement sensitivity matrices in which
%        each element is a matrix corresponding to a period in the forecast
%        horizon at time t = 1,2,...,numPeriods. If the lengths of the
%        observation vector y(t) and state vector x(t) are constant, then each
%        element is an n-by-m matrix; however, if the length of y(t) or x(t)
%        changes, then the elements are matrices of various sizes. The length
%        of the cell array must be at least numPeriods, and any elements
%        beyond the forecast horizon are ignored. By default, the last
%        coefficient of the input SSM model Mdl is used for all future periods.
%
%    'D' Cell vector of forecasted observation innovation matrices in which
%        each element is a matrix corresponding to a period in the forecast
%        horizon at time t = 1,2,...,numPeriods. If the lengths of the
%        observation vector y(t) and innovation vector e(t) are constant, then
%        each element is an n-by-h matrix; however, if the length of y(t) or
%        e(t) changes, then the elements are matrices of various sizes. The
%        length of the cell array must be at least numPeriods, and any
%        elements beyond the forecast horizon are ignored. By default, the
%        last coefficient of the input SSM model Mdl is used for all future
%        periods.
%
%    'Predictors0' T-by-d matrix of common predictor variables used to
%                  include a regression component in the observation equation.
%                  Observations at time t are deflated such that
%
%                  [y(t) - z(t)*b] = C * x(t) + D * e(t)
%
%                  where z(t) is a vector of predictor variables and b is
%                  the regression coefficient vector (see below). The default
%                  is an empty matrix (no regression component)
%
%    'PredictorsF' numPeriods-by-d matrix of common predictor variables used to
%                  include a regression component in the observation equation.
%
%    'Beta'        d-by-n matrix of regression coefficients associated with
%                  predictors (see above).
%
%  Output Arguments:
%
%    Y - Point forecasts of observations, E[y(t)|y(t-1),...,y(1)], for
%      t = 1,2,...,numPeriods. For time-invariant models in which the length
%      of each observation vector (n) is the same, this is a numPeriods-by-n
%      matrix. For time-varying models in which the length of the observation
%      vector changes, this is a numPeriods-by-1 cell array in which each
%      element contains a time-varying n-element vector of forecasts associated
%      with the corresponding period.
%
%    YMSE - Forecast error variances of future observations,
%      Cov[y(t)|y(t-1),...,y(1)], for t = 1,2,...,numPeriods. For
%      time-invariant models in which the length of each observation vector
%      (n) is the same, this is a numPeriods-by-n matrix. For time-varying
%      models in which the length of the observation vector changes, this is
%      a numPeriods-by-1 cell array in which each element contains a
%      time-varying n-element vector of forecast error variances associated
%      with the corresponding period.
%
%    X - Point forecasts of states, E[x(t)|y(t-1),...,y(1)], for t = 1,2,...,
%      numPeriods. For time-invariant models in which the length of each state
%      vector (m) is the same, this is a numPeriods-by-m matrix. For
%      time-varying models in which the length of the state vector changes,
%      this is a numPeriods-by-1 cell array in which each element contains a
%      time-varying m-element vector of forecasts associated with the
%      corresponding period.
%
%    XMSE - Forecast error variances of future states,
%      Cov[x(t)|y(t-1),...,y(1)], for t = 1,2,...,numPeriods. For
%      time-invariant models in which the length of each state vector (m) is
%      the same, this is a numPeriods-by-m matrix. For time-varying models
%      in which the length of the state vector changes, this is a
%      numPeriods-by-1 cell array in which each element contains a time-varying
%      m-element vector of forecast error variances associated with the
%      corresponding period.
%
%  See also SSM, FILTER, SMOOTH, SIMULATE, ESTIMATE.
%
%    Other functions named forecast
%
%       abstvar/forecast         garch/forecast
%       arima/forecast           generic/forecast
%       conjugateblm/forecast    gjr/forecast
%       customblm/forecast       regARIMA/forecast
%       diffuseblm/forecast      semiconjugateblm/forecast
%       egarch/forecast          varm/forecast
%       empiricalblm/forecast    vecm/forecast
%

author: jmaih


text:
%% Machine Learning Online Class - Exercise 3 | Part 1: One-vs-all

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     lrCostFunction.m (logistic regression cost function)
%     oneVsAll.m
%     predictOneVsAll.m
%     predict.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Setup the parameters you will use for this part of the exercise
input_layer_size  = 400;  % 20x20 Input Images of Digits
num_labels = 10;          % 10 labels, from 1 to 10
                          % (note that we have mapped "0" to label 10)

%% =========== Part 1: Loading and Visualizing Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  You will be working with a dataset that contains handwritten digits.
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

load('ex3data1.mat'); % training data stored in arrays X, y
m = size(X, 1);

% Randomly select 100 data points to display
rand_indices = randperm(m);
sel = X(rand_indices(1:100), :);

displayData(sel);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ============ Part 2a: Vectorize Logistic Regression ============
%  In this part of the exercise, you will reuse your logistic regression
%  code from the last exercise. You task here is to make sure that your
%  regularized logistic regression implementation is vectorized. After
%  that, you will implement one-vs-all classification for the handwritten
%  digit dataset.
%

% Test case for lrCostFunction
fprintf('\nTesting lrCostFunction() with regularization');

theta_t = [-2; -1; 1; 2];
X_t = [ones(5,1) reshape(1:15,5,3)/10];
y_t = ([1;0;1;0;1] >= 0.5);
lambda_t = 3;
[J, grad] = lrCostFunction(theta_t, X_t, y_t, lambda_t);

fprintf('\nCost: %f\n', J);
fprintf('Expected cost: 2.534819\n');
fprintf('Gradients:\n');
fprintf(' %f \n', grad);
fprintf('Expected gradients:\n');
fprintf(' 0.146561\n -0.548558\n 0.724722\n 1.398003\n');

fprintf('Program paused. Press enter to continue.\n');
pause;
%% ============ Part 2b: One-vs-All Training ============
fprintf('\nTraining One-vs-All Logistic Regression...\n');

lambda = 0.1;
[all_theta] = oneVsAll(X, y, num_labels, lambda);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 3: Predict for One-Vs-All ================

pred = predictOneVsAll(all_theta, X);

fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);

repo: machine-learning


text:
function [me,ASAcontrol] = moderr_e(varargin)
%MODERR_E Model error
%   ME = MODERR_E(AR_EST,MA_EST,AR_REF,MA_REF,N_OBS) is the model error
%   of an ARMA model estimated from N_OBS observations, with parameter
%   vectors AR_EST and MA_EST, applied for modeling a reference process
%   with ARMA parameters AR_REF and MA_REF. It is a measure for the
%   accuracy of the estimated model. The lower the model error, the
%   better is the model accuracy. For estimated models, its asymptotic
%   minimum value equals the AR-order plus the MA-order of the reference
%   process.
%
%   MODERR_E is an ARMASA main function.

%   References: P. M. T. Broersen, The Quality of Models for ARMA
%               Processes, IEEE Transactions on Signal Processing,
%               Vol. 46, No. 6,, June 1998, pp. 1749-1752.

%Header
%==============================================================================

%Declaration of variables
%------------------------

%Declare and assign values to local variables
%according to the input argument pattern
[ar_est,ma_est,ar_ref,ma_ref,n_obs,ASAcontrol]=ASAarg(varargin, ...
{'ar_est'   ;'ma_est'   ;'ar_ref'   ;'ma_ref'   ;'n_obs'    ;'ASAcontrol'}, ...
{'isnumeric';'isnumeric';'isnumeric';'isnumeric';'isnumeric';'isstruct'  }, ...
{'ar_est'   ;'ma_est'   ;'ar_ref'   ;'ma_ref'   ;'n_obs'                 });

if isequal(nargin,1) & ~isempty(ASAcontrol)
      %ASAcontrol is the only input argument
   ASAcontrol.error_chk = 0;
   ASAcontrol.run = 0;
end

%ARMASA-function version information
%-----------------------------------

%This ARMASA-function is characterized by
%its current version,
ASAcontrol.is_version = [2000 12 30 20 0 0];
%and its compatability with versions down to,
ASAcontrol.comp_version = [2000 12 30 20 0 0];

%This function calls other functions of the ARMASA
%toolbox. The versions of these other functions must
%be greater than or equal to:
ASAcontrol.req_version.convol_e = [2000 12 13 21 0 0];
ASAcontrol.req_version.arma2cor_e = [2000 12 30 20 0 0];

%Checks
%------

if ~isfield(ASAcontrol,'error_chk') | ASAcontrol.error_chk
      %Perform standard error checks
   %Input argument format checks
   ASAcontrol.error_chk = 1;
   if ~isnum(ar_est)
      error(ASAerr(11,'ar_est'))
   elseif ~isvector(ar_est)
      error(ASAerr(15,'ar_est'))
   elseif size(ar_est,1)>1
      ar_est = ar_est';
      warning(ASAwarn(25,{'column';'ar_est';'row'},ASAcontrol))
   end
   if ~isnum(ma_est)
      error(ASAerr(11,'ma_est'))
   elseif ~isvector(ma_est)
      error(ASAerr(15,'ma_est'))
   elseif size(ma_est,1)>1
      ma_est = ma_est';
      warning(ASAwarn(25,{'column';'ma_est';'row'},ASAcontrol))
   end
   if ~isnum(ar_ref)
      error(ASAerr(11,'ar_ref'))
   elseif ~isvector(ar_ref)
      error(ASAerr(15,'ar_ref'))
   elseif size(ar_ref,1)>1
      ar_ref = ar_ref';
      warning(ASAwarn(25,{'column';'ar_ref';'row'},ASAcontrol))
   end
   if ~isnum(ma_ref)
      error(ASAerr(11,'ma_ref'))
   elseif ~isvector(ma_ref)
      error(ASAerr(15,'ma_ref'))
   elseif size(ma_ref,1)>1
      ma_ref = ma_ref';
      warning(ASAwarn(25,{'column';'ma_ref';'row'},ASAcontrol))
   end
   if ~isnum(n_obs) | ~isintscalar(n_obs)
      error(ASAerr(17,'n_obs'))
   end

   %Input argument value checks
   if ~(isreal(ar_est) & isreal(ma_est) & ...
         isreal(ar_ref) & isreal(ma_ref))
      error(ASAerr(13))
   end
   if ar_est(1)~=1
      error(ASAerr(23,{'ar_est','parameter'}))
   end
   if ma_est(1)~=1
      error(ASAerr(23,{'ma_est','parameter'}))
   end
   if ar_ref(1)~=1
      error(ASAerr(23,{'ar_ref','parameter'}))
   end
   if ma_ref(1)~=1
      error(ASAerr(23,{'ma_ref','parameter'}))
   end
   if n_obs<1
      error(ASAerr(41,'n_obs'))
   end
end

if ~isfield(ASAcontrol,'version_chk') | ASAcontrol.version_chk
      %Perform version check
   ASAcontrol.version_chk = 1;

   %Make sure the requested version of this function
   %complies with its actual version
   ASAversionchk(ASAcontrol);

   %Make sure the requested versions of the called
   %functions comply with their actual versions
   arma2cor_e(ASAcontrol);
   convol_e(ASAcontrol);
end

if ~isfield(ASAcontrol,'run') | ASAcontrol.run
   ASAcontrol.run = 1;
end

if ASAcontrol.run %Run the computational kernel
   ASAcontrol.version_chk = 0;
   ASAcontrol.error_chk = 0;

%Main
%=====================================================

%Determine the model error by evaluating the power
%gain of a modified ARMA model
[cor,gain] = arma2cor_e...
   (convol_e(ar_ref,ma_est,ASAcontrol),...
    convol_e(ar_est,ma_ref,ASAcontrol),0,ASAcontrol);
me = n_obs*(gain-1);

%Footer
%=====================================================

else %Skip the computational kernel
   %Return ASAcontrol as the first output argument
   if nargout>1
      warning(ASAwarn(9,mfilename,ASAcontrol))
   end
   me = ASAcontrol;
   ASAcontrol = [];
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% former versions        P.M.T. Broersen        broersen@tn.tudelft.nl
% [2000 12 30 20 0 0]    W. Wunderink           wwunderink01@freeler.nl

author: Sable

repo: mcbench-benchmarks