text:
function eval(obj, inputs, derOutputs, stopLayer)
%EVAL Evaluate the DAGNN
%   EVAL(obj, inputs) evaluates the DaG for the specified input
%   values. `inputs` is a cell array of the type `{'inputName',
%   inputValue, ...}`. This call results in a forward pass through the
%   graph, computing the values of the output variables. These can
%   then be accessed using the `obj.vars(outputIndex)` property of the
%   DaG object. The index of an output can be obtained using the
%   `obj.getOutputIndex(outputName)` call.
%
%   EVAL(obj, inputs, derOutputs) evaluates the DaG forward and then
%   backward, performing backpropagation. Similar to `inputs`,
%   `derOutputs` is a cell array of the type {'outputName',
%   outputDerValue, ...} of output derivatives.
%
%   # Understanding backpropagation
%
%   Only those outputs for which an `outputDerValue` which is
%   non-empty are involved in backpropagation, while the others are
%   ignored. This is useful to attach to the graph auxiliary layers to
%   compute errors or other statistics, without however involving them
%   in backpropagation.
%
%   Usually one starts backpropagation from scalar outptus,
%   corresponding to loss functions. In this case `outputDerValue` can
%   be interpreted as the weight of that output and is usually set to
%   one. For example: `{'objective', 1}` backpropagates from the
%   `'objective'` output variable with a weight of 1.
%
%   However, in some cases the DaG may contain more than one such
%   node, for example because one has more than one loss function.  In
%   this case `{'objective1', w1, 'objective2', w2, ...}` allows to
%   balance the different objectives.
%
%   Finally, one can backpropagate from outputs that are *not*
%   scalars. While this is unusual, it is possible by specifying a
%   value of `outputDerValue` that has the same dimensionality as the
%   output; in this case, this value is used as a matrix of weights,
%   or projection.
%
%   # Factors affecting evaluation
%
%   There are several factors affecting evaluation:
%
%   * The *evaluation mode* can be either `normal` or `test`. Layers
%     may behave differently depending on the mode. For example,
%     dropout becomes a pass-through layer in test mode and batch
%     normalization use fixed moments (this usually improves the test
%     performance significantly).
%
%   * By default, the DaG aggressively conserves memory. This is
%     particularly important on the GPU, where memory is
%     scarce. However, this also means that the values of most
%     variables and of their derivatives are dropped during the
%     computation. For debugging purposes, it may be interesting to
%     observe these variables; in this case you can set the
%     `obj.conserveMemory` property of the DaG to `false`. It is also
%     possible to preserve individual variables by setting the
%     property `obj.vars(v).precious` to `true`.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

obj.computingDerivative = nargin > 2 && ~isempty(derOutputs) ;

if ~iscell(inputs), error('INPUTS is not a cell array.') ; end
if obj.computingDerivative && ~iscell(derOutputs), error('DEROUTPUTS is not a cell array.') ; end

if nargin < 4
  stopLayer = [];
end
% -------------------------------------------------------------------------
% Forward pass
% -------------------------------------------------------------------------

% set the input values
v = obj.getVarIndex(inputs(1:2:end)) ;
if any(isnan(v))
  broken = find(isnan(v)) ;
%   error('No variable of name ''%s'' could be found in the DAG.', inputs{2*broken(1)-1}) ;
end

for i = 1:2:numel(inputs)
  v = obj.getVarIndex(inputs{i}) ;
  if ~isnan(v)
      switch obj.device
        case 'cpu', obj.vars(v).value = gather(inputs{i+1}) ;
        case 'gpu', obj.vars(v).value = gpuArray(inputs{i+1}) ;
      end
  end
end

inputs = [] ;
obj.numPendingVarRefs = [obj.vars.fanout] ;
for l = obj.executionOrder
  time = tic ;
  obj.layers(l).block.forwardAdvanced(obj.layers(l)) ;
  obj.layers(l).forwardTime = toc(time) ;

  if stopLayer == l
    break;
  end
end

% -------------------------------------------------------------------------
% Backward pass
% -------------------------------------------------------------------------

if ~obj.computingDerivative, return ; end

% conserve memory
if obj.backpropDepth
for l = obj.executionOrder
  [obj.vars(obj.layers(l).inputIndexes).value] = deal([]) ;
  if strcmp(obj.layers(l).name, obj.backpropDepth)
    break;
  end
end
end

v = obj.getVarIndex(derOutputs(1:2:end)) ;
[obj.vars(v).der] = deal(derOutputs{2:2:end}) ;
derOutputs = [] ;

obj.numPendingVarRefs = zeros(1, numel(obj.vars)) ;
obj.numPendingParamRefs = zeros(1, numel(obj.params)) ;
for l = fliplr(obj.executionOrder)
  time = tic ;
  if strcmp(obj.layers(l).name, obj.backpropDepth)
    [obj.vars.value] = deal([]) ;
    [obj.vars.der] = deal([]) ;
    break;
  end
  obj.layers(l).block.backwardAdvanced(obj.layers(l)) ;
  obj.layers(l).backwardTime = toc(time) ;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function pdf = dipole_pdf ( x, a, b )

%*****************************************************************************80
%
%% DIPOLE_PDF evaluates the Dipole PDF.
%
%  Discussion:
%
%    PDF(X)(A,B) =
%        1 / ( PI * ( 1 + X**2 ) )
%      + B**2 * ( ( 1 - X**2 ) * cos ( 2 * A ) + 2.0D+00 * X * sin ( 2 * A ) )
%      / ( PI * ( 1 + X )**2 )
%
%    Densities of this kind commonly occur in the analysis of resonant
%    scattering of elementary particles.
%
%    DIPOLE_PDF(X)(A,0) = CAUCHY_PDF(X)(A)
%    A = 0, B = 1 yields the single channel dipole distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Knop,
%    Algorithm 441,
%    ACM Transactions on Mathematical Software.
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    A is arbitrary, but represents an angle, so only 0 <= A <= 2 * PI
%      is interesting,
%    and -1.0D+00 <= B <= 1.0.
%
%    Output, real PDF, the value of the PDF.
%
  pdf = 1.0 / ( pi * ( 1.0 + x * x ) ) ...
    + b * b * ( ( 1.0 - x * x ) * cos ( 2.0 * a ) ...
    + 2.0 * x * sin ( 2.0 * x ) ) / ( pi * ( 1.0 + x * x ) * ( 1.0 + x * x ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%

sim = exp(-sum((x1-x2).^2)/2/(sigma^2));

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

end

author: zzlyw

repo: machine-learning-exercises


text:
function r= cauchyrnd(varargin)

% USAGE:       r= cauchyrnd(a, b, n, ...)
%
% Generate random numbers from the Cauchy distribution, r= a + b*tan(pi*(rand(n)-0.5)).
%
% ARGUMENTS:
% a (default value: 0.0) must be scalars or size(x).
% b (b>0, default value: 1.0) must be scalars or size(x).
% n and onwards (default value: 1) specifies the dimension of the output.
%
% EXAMPLE:
% r= cauchyrnd(0, 1, 10); % A 10 by 10 array of random values, Cauchy distributed.
%
% SEE ALSO:    cauchycdf, cauchyfit, cauchyinv, cauchypdf.
%
% Copyright (C) Peder Axensten <peder at axensten dot se>
%
% HISTORY:
% Version 1.0, 2006-07-10.
% Version 1.1, 2006-07-26.
% - Added cauchyfit to the cauchy package.
% Version 1.2, 2006-07-31:
% - cauchyinv(0, ...) returned a large negative number but should be -Inf.
% - Size comparison in argument check didn't work.
% - Various other improvements to check list.

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

	% Default values
	a=	0.0;
	b=	1.0;
	n=	1;

	% Check the arguments
	if(nargin >= 1)
		a=	varargin{1};
		if(nargin >= 2)
			b=			varargin{2};
			b(b <= 0)=	NaN;	% Make NaN of out of range values.
			if(nargin >= 3),	n=	[varargin{3:end}];		end
		end
	end

	% Generate
	r=	cauchyinv(rand(n), a, b);
end

author: Sable

repo: mcbench-benchmarks


text:
function enabled = isEnabled(package)
%  enabled = isEnabled(package)
%
% Takes
%
% 1. package - a string representing an SBML L3 package
%
% Returns
%
% 1. enabled
%   - 1, if the package can be used by this instance of libSBML
%   - 0, otherwise

%<!---------------------------------------------------------------------------
% This file is part of libSBML.  Please visit http://sbml.org for more
% information about SBML, and the latest version of libSBML.
%
% Copyright (C) 2013-2018 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%     3. University of Heidelberg, Heidelberg, Germany
%
% Copyright (C) 2009-2013 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%
% Copyright (C) 2006-2008 by the California Institute of Technology,
%     Pasadena, CA, USA
%
% Copyright (C) 2002-2005 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. Japan Science and Technology Agency, Japan
%
% This library is free software; you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation.  A copy of the license
% agreement is provided in the file named "LICENSE.txt" included with
% this software distribution and also available online as
% http://sbml.org/software/libsbml/license.html
%----------------------------------------------------------------------- -->

% assume not enabled
enabled = 0;

supported = {'fbc', 'qual', 'groups'};
if ~ischar(package)
    disp('argument must be a string representing an SBML L3 package');
elseif ~ismember(supported, package)
    return;
end;

if (isoctave() == '0')
  filename = fullfile(tempdir, 'pkg.xml');
else
  filename = fullfile(pwd, 'pkg.xml');
end;
writeTempFile(filename, package);

version_field = strcat(package, '_version');
try
  [m, e] = TranslateSBML(filename, 1, 0);

  if (isempty(e) && isfield(m, version_field) == 1 )
    enabled = 1;
  end;

  delete(filename);

catch err
  disp(err.identifier)
  delete(filename);

  return
end;

function writeTempFile(filename, package)

fout = fopen(filename, 'w');

fprintf(fout, '<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n');
fprintf(fout, '<sbml xmlns=\"http://www.sbml.org/sbml/level3/version1/core\" ');
fprintf(fout, strcat('xmlns:', package, '=\"http://www.sbml.org/sbml/level3/version1/', package, '/version1\"\n'));
fprintf(fout, 'level=\"3\" version=\"1\"\n');
fprintf(fout, strcat(package, ':required=\"false\">\n'));
fprintf(fout, '  <model/>\n</sbml>\n');

fclose(fout);

author: opencobra

repo: cobratoolbox


text:
function poly_cof = r8poly_basis ( ntab, xtab )

%*****************************************************************************80
%
%% R8POLY_BASIS computes all Lagrange basis polynomial in standard form.
%
%  Discussion:
%
%    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
%    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
%    XTAB(J) for J not equal to I, and 1 when J is equal to I.
%
%    The Lagrange basis polynomials have the property that the interpolating
%    polynomial through a set of NTAB data points (XTAB,YTAB) may be
%    represented as
%
%      P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X)
%
%    Higher order interpolation at selected points may be accomplished
%    using repeated X values, and scaled derivative values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NTAB, the number of data points XTAB.
%
%    Input, real XTAB(NTAB), the X values upon which the
%    Lagrange basis polynomial is to be based.
%
%    Output, real POLY_COF(NTAB,NTAB), the polynomial
%    coefficients for the I-th Lagrange basis polynomial are stored
%    in column I.  POLY_COF(1,I) is the constant term, and POLY_COF(1,NTAB)
%    is the coefficient of X**(NTAB-1).
%

%
%  Initialize POLY_COF to the identity matrix.
%
  poly_cof(1:ntab,1:ntab) = 0.0;
  for i = 1 : ntab
    poly_cof(i,i) = 1.0;
  end
%
%  Compute the divided difference table for the IVAL-th Lagrange basis
%  polynomial.
%
  for i = 1 : ntab
    poly_cof(1:ntab,i) = ( data_to_dif ( ntab, xtab, poly_cof(1:ntab,i) ) )';
  end
%
%  Convert the divided difference table coefficients to standard polynomial
%  coefficients.
%
  for i = 1 : ntab
    poly_cof(1:ntab,i) = ( dif_to_r8poly ( ntab, xtab, poly_cof(1:ntab,i) ) )';
  end

  return
end


text:
 function ob = ctranspose(ob)
%function ob = ctranspose(ob)
%	"ctranspose" method for Gtomo2 class

%	transpose base object
base = ob.base;
ob.base = base';

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [V,cc]=intstab(G,Gc,H)
G1=zpk(minreal(feedback(Gc*G,H)));
G2=zpk(Gc*G*H);
V=0; cc=[];
vv=G1.p{1}; ii=find(real(vv)>=0);
if length(ii)>0, V=1; cc=vv(ii);
else
   z=G2.z{1}; p=G2.p{1};
   for i=1:length(z),
      ii=find(abs(z(i)-p)<1000*eps);
      if length(ii)>0
         if real(z(i))>=0,
            V=2; cc=[cc,z(i)]; break;
end, end, end, end


text:
%==================================================================
%                        General Data File
% Title: TETRAHEDRA
% Units: SI
% Dimensions: 2D
% Type of problem: Plane_Stress
% Type of Phisics: ELASTIC
% Micro/Macro: MACRO
%
%==================================================================

%% Data

Data_prb = {
'TRIANGLE';
'SI';
'2D';
'Plane_Stress';
'THERMAL';
'MACRO';
};

%% Coordinates
% Node                X                Y                Z

coord = [
   1        00000        00000
     2  1.00000e-02  1.00000e-02
     3        00000  2.00000e-02
     4  2.00000e-02        00000
     5  2.00000e-02  2.00000e-02
     6  1.00000e-02  3.00000e-02
     7  3.00000e-02  1.00000e-02
     8        00000  4.00000e-02
     9  4.00000e-02        00000
    10  3.00000e-02  3.00000e-02
    11  4.00000e-02  2.00000e-02
    12  2.00000e-02  4.00000e-02
    13  4.00000e-02  4.00000e-02
];

%% Conectivities
% Element        Node(1)                Node(2)                Node(3)                Node(4)                Material

connec = [
     1      5      2      4
     2      3      2      5
     3      1      2      3
     4      4      2      1
     5     12      6      5
     6      8      6     12
     7      3      6      8
     8      5      6      3
     9     11      7      9
    10      5      7     11
    11      4      7      5
    12      9      7      4
    13     13     10     11
    14     12     10     13
    15      5     10     12
    16     11     10      5
];

%% Variable Prescribed
% Node            Dimension                Value

dirichlet_data = [
                            1 1 0
                            3 1 0
                            8 1 0
];

%% Force Prescribed
% Node                Dimension                Value

pointload_complete = [
11 1 -1
];

%% Volumetric Force
% Element        Dim                Force_Dim

Vol_force = [
];

%% Group Elements
% Element        Group_num

Group = [
];

%% Initial Holes
% Elements that are considered holes initially
% Element

Initial_holes = [
];

%% Boundary Elements
% Elements that can not be removed
% Element

Boundary_elements = [
];

%% Micro gauss post
%
% Element

Micro_gauss_post = [
];

%% Micro Slave-Master
% Nodes that are Slaves
% Nodes             Value (1-Slave,0-Master)

Micro_slave = [
];

%% Nodes solid
% Nodes that must remain
% Nodes

nodesolid = unique(pointload_complete(:,1));

%% External border Elements
% Detect the elements that define the edge of the domain
% Element               Node(1)           Node(2)

External_border_elements = [
];

%% External border Nodes
% Detect the nodes that define the edge of the domain
% Node

External_border_nodes = [
];

%% Materials
% Materials that have been used
% Material_Num              Mat_density        Young_Modulus        Poisson

Materials = [
];

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function Ath = threshold_components(A,options)

% post processing of spatial components
% for each component perform the following:
%   (i)     perform median filtering
%   (ii)    threshold components
%   (iii)   perform morphological closing
%   (iv)    extract largest connected component

% The thresholding can be determined in two ways:
% options.thr_method = 'max': For every component pixel values below
%       options.thr_method are discarded
% options.thr_method = 'nrg': The boundary is drawn around the value above
%        which a specified fraction of energy is explained (default 99%)

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

    defoptions = CNMFSetParms;
    if nargin < 2 || isempty(options); options = defoptions; end
%     defoptions.nrgthr = 0.9999;              % energy threshold
%     defoptions.clos_op = strel('square',3);  % morphological operator for closing
%     defoptions.medw = [3,3];                 % size of median filter

    if ~isfield(options,'thr_method') || isempty(options.thr_method); options.thr_method = defoptions.thr_method; end
    if ~isfield(options,'nrgthr') || isempty(options.nrgthr); options.nrgthr = defoptions.nrgthr; end
    if ~isfield(options,'maxthr') || isempty(options.maxthr); options.maxthr = defoptions.maxthr; end
    if ~isfield(options,'clos_op') || isempty(options.clos_op); options.clos_op = defoptions.clos_op; end
    if ~isfield(options,'conn_comp') || isempty(options.conn_comp); options.conn_comp = defoptions.conn_comp; end % extract largest connected component
    if ~isfield(options,'medw') || isempty(options.medw); options.medw = defoptions.medw; end
    if ~isfield(options,'d3') || isempty(options.d3); options.d3 = 1; end

    if options.d3 > 1; options.thr_method = 'nrg'; end
    [d,nr] = size(A);
    Ath = spalloc(d,nr,nnz(A));
    Ath(:,nr-options.nb+1:nr) = A(:,nr-options.nb+1:nr);
    indf = cell(nr,1);
    valf = cell(nr,1);
    parfor i = 1:nr-options.nb
        A_temp = reshape(full(A(:,i)),options.d1,options.d2,options.d3);
        for z = 1:options.d3
            A_temp(:,:,z) = medfilt2(A_temp(:,:,z),options.medw);
        end
        if strcmpi(options.thr_method,'nrg')
            A_temp = A_temp(:);
            [temp,ind] = sort(A_temp(:).^2,'ascend');
            temp =  cumsum(temp);
            ff = find(temp > (1-options.nrgthr)*temp(end),1,'first');
            BW = zeros(options.d1,options.d2,options.d3);
            BW(ind(ff:d)) = 1;
            for z = 1:options.d3
                BW(:,:,z) = imclose(BW(:,:,z),options.clos_op);
            end
            [L,NUM] = bwlabeln(BW,8*(options.d3==1) + 6*(options.d3~=1));
            if NUM > 0
                nrg = zeros(NUM,1);
                for l = 1:NUM
                    ff = (L==l);
                    nrg(l) = sum(A_temp(ff).^2);
                end
                [~,indm] = max(nrg);
                ff = find(L==indm);
                indf{i} = ff;
                valf{i} = A_temp(ff);
            else
                valf{i} = 0;
            end
        else
            if ~strcmpi(options.thr_method,'max')
                fprintf('Unknown thresholding method, using thresholding based on max value. \n');
            end
            A_temp(A_temp<options.maxthr*max(A_temp(:))) = 0;
            BW = imclose(A_temp>0,options.clos_op);
            if options.conn_comp;
                BW = bwareafilt(BW,1);
            end
            ff = find(BW>0);
            indf{i} = ff;
            if ~isempty(ff)
                valf{i} = A_temp(ff);
            else
                valf{i} = 0;
            end
        end
    end
    for i = 1:nr-options.nb
        if ~isempty(indf{i})
            Ath(indf{i},i) = valf{i};
        end
    end
end

author: flatironinstitute

repo: CaImAn-MATLAB


function [ L, D, C, Cinv ] = get_transform_matrices_2d( R, method )
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here
if method == 1
    % LDL Without the permutation
    [L, D] = ldl(R);
    C = inv(L);
    Cinv = L;

elseif method == 2
    % Eigen Decomp
    [V,D] = eig(R);
    C = V.';
    Cinv = V;

elseif method == 3
    % Cholesky Decomp (Special case of LDL, D=I)
    L = chol(R);  % NOTE: this is the UPPER triangular transform
    L = L.';
    C = inv(L);
    Cinv = L;
    D = eye(2,2);

elseif method == 4  %%% ANALYTICAL
    rho = R(1,2);
    L = [1 0; rho 1];
    D = diag([1, 1 - rho^2]);
    C = [1 0; -rho 1];
    % C = inv(L);
    Cinv = L;
end

end


text:
%compute yinflection
function [data,units]=compute_yinflection_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
yinflection_mm=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    sz = size(trx(larva).yspine_mm);
    yinflection_mm{1,i} = trx(larva).yspine_mm(sub2ind(sz,trx(larva).inflectionpointdistance,1:sz(2)));
%
%     yinflection_mm{1,i}=zeros(1,size(trx(larva).yspine_mm,2));
%     for j=1:size(trx(larva).yspine_mm,2)
%     yinflection_mm{1,i}(j)=trx(larva).yspine_mm(trx(larva).inflectionpointdistance(j),j);
%     end
end

units=parseunits('mm');
data=yinflection_mm;

author: kristinbranson

repo: JAABA


text:
%% Example 8.6: Simulating from a trigonometric nonlinear SDE
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%% Gaussian approximation

  % Lock random seed
  if exist('rng') % Octave doesn't have rng
    rng(1,'twister')
  else
    randn('state',1);
    rand('state',1);
  end

  % Parameters
  tspan = 0:1:10;

  % The model
  f = @(x,t) -(1/10)^2*sin(x).*cos(x).^3;
  L = @(x,t) 1/10*cos(x).^2;

  % The derivatives
  df  = @(x,t) -1/100*cos(x).^2.*(2*cos(2*x)-1);
  ddf = @(x,t) 1/100*(sin(2*x) + 2*sin(4*x));
  dL  = @(x,t) -1/5*sin(x).*cos(x);
  ddL = @(x,t) -1/5*cos(2*x);

  % Exact solution
  solfun = @(wt,x0) atan(1/10*wt +tan(x0));

  % Intial
  x0 = 1;

%% Sample

  x = zeros(1,10000);
  xem   = x;
  xw20  = x;
  xw20g = x;

  % Lock random seed
  if exist('rng') % Octave doesn't have rng
    rng(1,'twister')
  else
    randn('state',1);
    rand('state',1);
  end

  for j=1:size(x,2)

    % Use Euler-Maruyama
    foo = eulermaruyama_weak(f,L,tspan,x0,1);
    xem(:,j) = foo(:,end);

    % Report
    if rem(j,100)==0, j, end

  end

  % Lock random seed
  if exist('rng') % Octave doesn't have rng
    rng(1,'twister')
  else
    randn('state',1);
    rand('state',1);
  end

  for j=1:size(x,2)

    % Weak order 2.0
    foo = w20scalar({f,df,ddf},{L,dL,ddL},tspan,x0,1,false);
    xw20(:,j) = foo(:,end);

    % Report
    if rem(j,100)==0, j, end

  end

  % Lock random seed
  if exist('rng') % Octave doesn't have rng
    rng(1,'twister')
  else
    randn('state',1);
    rand('state',1);
  end

  for j=1:size(x,2)

    % Weak order 2.0 (Gaussian increments)
    foo = w20scalar({f,df,ddf},{L,dL,ddL},tspan,x0,1,true);
    xw20g(:,j) = foo(:,end);

    % Store
    x(:,j) = foo(:,end);

    % Report
    if rem(j,100)==0, j, end

  end

%% Visualize

  % Lock random seed
  if exist('rng') % Octave doesn't have rng
    rng(1,'twister')
  else
    randn('state',1);
    rand('state',1);
  end

  % Samples from the exact solution
  xe = solfun(sqrt(tspan(end))*randn(1,200000),x0);

  % Bins
  nbins = 64;
  t = linspace(min(xe),max(xe)+.1,nbins);

  figure(2); clf; hold on

    % Show solution
    n = histc(xe,t);
    fill(t,n/numel(xe),1,'FaceColor',[.7 .7 .7],'EdgeColor',[.7 .7 .7])

    % Show solution w2.0
    n = histc(xw20,t);
    stairs(t-(t(2)-t(1))/2,n/numel(xw20),'-k')

    % Limits
    xlim([.35 1.25])
    lims = ylim;

    % Label
    xlabel('$x$')

    % Ticks
    set(gca,'XTick',0:.2:1.2)

  figure(3); clf; hold on

    % Show solution
    n = histc(xe,t);
    fill(t,n/numel(xe),1,'FaceColor',[.7 .7 .7],'EdgeColor',[.7 .7 .7])

    % Show solution w2.0 (Gaussian increments)
    n = histc(xw20g,t);
    stairs(t-(t(2)-t(1))/2,n/numel(xw20g),'-k')

    % Limits
    %xlim([min(t) max(t)])
    xlim([.35 1.25])
    ylim(lims)

    % Show legend
    legend('Exact','Weak order $2.0$')

    % Label
    xlabel('$x$')

    % Ticks
    set(gca,'XTick',0:.2:1.2)

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
% Jiao Xianjun (putaoshu@msn.com; putaoshu@gmail.com)
% analyze PSS detection under big frequency offset
% A script of project: https://github.com/JiaoXianjun/rtl-sdr-LTE

% according to simulation result, gold configuration is that use whole PSS (length 137), and fo step 50kHz
% simulation result:
% len pss 137 fo step 50kHz SNR 0dB 10000 cost 29.9569s fail rate 0.003
% len pss 137 fo step 50kHz SNR 3dB 10000 cost 27.4461s fail rate 0
% len pss 137 fo step 25kHz SNR 0dB 10000 cost 40.6559s fail rate 0.0025
% len pss 137 fo step 25kHz SNR 3dB 10000 cost 41.6056s fail rate 0
% len pss 128 fo step 50kHz SNR 0dB 10000 cost 23.6994s fail rate 0.0091
% len pss 128 fo step 50kHz SNR 3dB 10000 cost 23.0074s fail rate 0.0002
% len pss 128 fo step 25kHz SNR 0dB 10000 cost 41.1991s fail rate 0.0078
% len pss 128 fo step 25kHz SNR 3dB 10000 cost 41.4224s fail rate 0

% function test_pss_fo_corr
clear all;
close all;

[fd_pss, td_pss] = pss_gen;

pss_idx = 1;
% pss = td_pss(end-127:end, pss_idx);
pss = td_pss(:, pss_idx);
len_pss = length(pss);

fo_search_set = [-200e3:50e3:200e3];
% fo_search_set = 0;
pss_set = kron(ones(1, length(fo_search_set)), pss).*exp(1i.*2.*pi.*(1./1.92e6).*(0:(length(pss)-1)).'*fo_search_set);
fo_step = (fo_search_set(2) - fo_search_set(1))/1e3;

var_pss = mean(abs(td_pss(:, pss_idx)).^2);

r = zeros(1, 2*137);
tmp_total = zeros(1, length(fo_search_set));
num_test = 10000;
num_fail = 0;
snr = 3;
sigma2 = var_pss/(10^(snr/10));
correct_pos = 147 - (length(pss)-128);
tic;
for idx = 1:num_test
    r_pss = [(randn(137,1)+1i.*randn(137,1)).*sqrt(var_pss./2); td_pss(:, pss_idx); (randn(137,1)+1i.*randn(137,1)).*sqrt(var_pss./2)];
    r_pss = r_pss + sqrt(sigma2/2).*(randn(length(r_pss),1)+1i.*randn(length(r_pss),1));

    fo = (2*rand-1)*200e3;
%     fo = (-7.5e3/1) + 0*15e3;
    r_pss_fo = r_pss.*exp(1i.*2.*pi.*fo.*(1./1.92e6).*(0:(length(r_pss)-1))');

    for i=1:(2*137)
        s = r_pss_fo(i:(i+length(pss)-1));

%         %  %---------method 1-----------------------
%         for j=1:length(fo_search_set)
%             tmp = conj(pss_set(:,j)).*s;
%             tmp = vec2mat(tmp, 8);
%             tmp = sum( abs(sum(tmp, 2)).^2 );
% %             r(i) = r(i) + tmp;
%             tmp_total(j) = tmp;
%         end
%         r(i) = max(tmp_total);
%         %  %---------end of method 1-----------------------

%         %  %---------method 2-----------------------
%         r(i) = abs(s(1:64)'*s(128:-1:65))^2;
%         %  %---------end of method 2-----------------------

        %  %---------method 3-----------------------
        tmp = conj(pss).*s;
        r(i) = max(abs(fft(tmp)).^2);
        %  %---------end of method 3-----------------------
    end
    plot(r); drawnow;
    % plot(138, r(138), 'rs');

    if mod(idx, 100) == 0
        disp(['len pss ' num2str(len_pss) ' fo step ' num2str(fo_step) 'kHz SNR ' num2str(snr) 'dB ' num2str(idx) ' cost ' num2str(toc) 's fail rate ' num2str(num_fail/idx) ]);
        tic;
    end

    [~, max_idx] = max(r);
    if max_idx ~= correct_pos
%         plot(r); hold on;
%         plot(correct_pos, r(correct_pos), 'rs');
        num_fail = num_fail + 1;
        disp(['idx ' num2str(idx) ' max_idx ' num2str(max_idx) ' fail rate ' num2str(num_fail/idx)]);
%         break;
    end
end

disp(['fail rate ' num2str(num_fail/idx)]);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function centroid = voronoi_polygon_centroid ( node, neighbor_num, ...
  neighbor_index, node_num, node_xy )

%*****************************************************************************80
%
%% VORONOI_POLYGON_CENTROID computes the centroid of a Voronoi polygon.
%
%  Discussion:
%
%    It is assumed that the Voronoi polygon is finite!  Every Voronoi
%    diagram includes some regions which are infinite, and for those,
%    this formula is not appropriate.
%
%    The routine is given the indices of the nodes that are neighbors of a
%    given "center" node.  A node is a neighbor of the center node if the
%    Voronoi polygons of the two nodes share an edge.  The triangles of the
%    Delaunay triangulation are formed from successive pairs of these neighbor
%    nodes along with the center node.
%
%    The assumption that the polygon is a Voronoi polygon is
%    used to determine the location of the boundaries of the polygon,
%    which are the perpendicular bisectors of the lines connecting
%    the center point to each of its neighbors.
%
%    The finiteness assumption is employed in part in the
%    assumption that the polygon is bounded by the finite
%    line segments from point 1 to 2, 2 to 3, ...,
%    M-1 to M, and M to 1, where M is the number of neighbors.
%
%    It is assumed that this subroutine is being called by a
%    process which has computed the Voronoi diagram of a large
%    set of nodes, so the arrays X and Y are dimensioned by
%    NODE_NUM, which may be much greater than the number of neighbor
%    nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, Sung Nok Chiu,
%    Spatial Tessellations: Concepts and Applications of Voronoi Diagrams,
%    Second Edition,
%    Wiley, 2000, page 490.
%
%  Parameters:
%
%    Input, integer NODE, the index of the node whose Voronoi
%    polygon is to be analyzed.  1 <= NODE <= NODE_NUM.
%
%    Input, integer NEIGHBOR_NUM, the number of neighbor nodes of
%    the given node.
%
%    Input, integer NEIGHBOR_INDEX(NEIGHBOR_NUM), the indices
%    of the neighbor nodes.  These indices are used to access the
%    X and Y arrays.  The neighbor nodes should be listed in the
%    (counter-clockwise) order in which they occur as one circles
%    the center node.
%
%    Input, integer NODE_NUM, the total number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Output, real CENTROID(2), the coordinates of the centroid
%    of the Voronoi polygon of node NODE.
%
  dim_num = 2;

  centroid(1:dim_num,1) = 0.0;

  if ( node < 1 | node_num < node )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VORONOI_POLYGON_CENTROID - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of input parameter NODE.\n' );
    error ( 'VORONOI_POLYGON_CENTROID - Fatal error!' )
  end

  pc(1:dim_num) = node_xy(1:dim_num,node)';

  i = neighbor_num;
  i = neighbor_index(i);

  pi(1:dim_num) = node_xy(1:dim_num,i)';

  j = 1;
  j = neighbor_index(j);

  pj(1:dim_num) = node_xy(1:dim_num,j)';

  a = ( pi(1) * pi(1) + pi(2) * pi(2) - pc(1) * pc(1) - pc(2) * pc(2) );
  b = ( pj(1) * pj(1) + pj(2) * pj(2) - pc(1) * pc(1) - pc(2) * pc(2) );
  c = 2.0 * ( ( pi(1) - pc(1) ) * ( pj(2) - pc(2) ) ...
            - ( pj(1) - pc(1) ) * ( pi(2) - pc(2) ) );
  uj(1) = ( a * ( pj(2) - pc(2) ) - b * ( pi(2) - pc(2) )  ) / c;
  uj(2) = ( a * ( pj(1) - pc(1) ) - b * ( pi(1) - pc(1) )  ) / c;

  for i = 1 : neighbor_num

    pi(1:dim_num) = pj(1:dim_num);
    ui(1:dim_num) = uj(1:dim_num);

    j = i4_wrap ( i+1, 1, neighbor_num );
    pj(1:dim_num) = node_xy(1:dim_num,j)';

    a = ( pi(1) * pi(1) + pi(2) * pi(2) - pc(1) * pc(1) - pc(2) * pc(2) );
    b = ( pj(1) * pj(1) + pj(2) * pj(2) - pc(1) * pc(1) - pc(2) * pc(2) );
    c = 2.0 * ( ( pi(1) - pc(1) ) * ( pj(2) - pc(2) ) ...
              - ( pj(1) - pc(1) ) * ( pi(2) - pc(2) ) );
    uj(1) = ( a * ( pj(2) - pc(2) ) - b * ( pi(2) - pc(2) )  ) / c;
    uj(2) = ( a * ( pj(1) - pc(1) ) - b * ( pi(1) - pc(1) )  ) / c;

    centroid(1) = centroid(1) + ( ui(2) - uj(2) ) ...
      * ( ( uj(1) + ui(1) )^2 - uj(1) * ui(1) );
    centroid(2) = centroid(2) + ( ui(1) - uj(1) ) ...
      * ( ( uj(2) + ui(2) )^2 - uj(2) * ui(2) );

  end

  area = voronoi_polygon_area ( node, neighbor_num, neighbor_index, ...
    node_num, node_xy );

  centroid(1:dim_num) = centroid(1:dim_num) / ( 6.0 * area );

  return
end


text:
function[h]=patchcontourf(varargin)
%PATCHCONTOURF  Generate filled contours using patches, with specified colors.
%
%   Frequently one wishes to make a filled contour plot of one field over
%   another with two different colormap schemes.  This is not currently
%   possible with CONTOURF because Matlab uses one colormap per axes.
%
%   For example, one may wish to plot sea surface temperature filled in
%   color with the continents filled in black.
%
%   PATCHCONTOURS gets around this problem by geneating filled contours
%   plots as patch objects, with out reference to the current colormap.
%   Colors of filled contours may be specified explicitly.
%
%   PATCHCONTOURF(Z,V) makes a filled contour plot of field Z at value V.
%   The interiors of the contour are colored black by default.
%
%   PATCHCONTOURF(Z,V,C) alternately uses color C as both the fill color
%   and the edge color.  C may be either a string, as in 'g' for green,
%   or a color array of length 3, such as [1/2 1/2 1/2] for gray.
%
%   PATCHCONTOURF(X,Y,Z,V) and PATCHCONTOURF(X,Y,Z,V,C) use vectors X and Y
%   to specify the contour locations.  X has length SIZE(X,2) while Y has
%   length SIZE(X,1).
%
%   H=PATCHCONTOURF(...) returns a handle to the patch object.
%
%   PATCHCONTOURF(...,'m_map') alternately works with R. Pawlowicz's M_MAP
%   mapping package, available at http://www.eos.ubc.ca/~rich/map.html.
%
%   'patchcontourf --f' generates a sample figure.
%
%   Usage: patchcontourf(z,v);
%          patchcontourf(z,v,c);
%          patchcontourf(x,y,z,v,c);
%          h=patchcontourf(x,y,z,v,c);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmp(varargin{1}, '--f')
    type makefigs_patchcontourf
    makefigs_patchcontourf;
    return
end

%Sort out string input arguments
str='k';
flag='mat';
if isstr(varargin{end})
    if length(varargin{end})>=3
        if strcmpi(varargin{end}(1:3),'m_m')||strcmpi(varargin{end}(1:3),'mat')
            flag=varargin{end};
            varargin=varargin(1:end-1);
        end
    end
end

if isstr(varargin{end})||length(varargin{end})==3
    str=varargin{end};
    varargin=varargin(1:end-1);
end

v=varargin{end};
varargin=varargin(1:end-1);

if length(varargin)==3
    x=varargin{1};
    y=varargin{2};
    z=varargin{3};
else
    x=[];
    y=[];
    z=varargin{1};
end

z1=zeros(size(z,1)+2,size(z,2)+2)+minmin(z);
z1(2:end-1,2:end-1)=z;

[ii,jj]=closedcurves(z1,v);

bool=ishold;
hold on

h=[];
if isempty(x)
    for i=1:length(ii)
        if strcmpi(flag(1:3),'m_m')
            m_patch(ii{i},jj{i},str,'edgecolor',str);
        else
            patch(ii{i},jj{i},str,'edgecolor',str);
        end
    end
else
    dx1=x(2)-x(1);
    dxN=x(end-1)-x(end);
    x=[x(1)-dx1;x(:);x(end)-dxN];

    dy1=y(2)-y(1);
    dyN=y(end-1)-y(end);
    y=[y(1)-dy1;y;y(end)-dyN];

    for i=1:length(ii)
        x1=interp1(1:length(x),x,ii{i},'linear');
        y1=interp1(1:length(y),y,jj{i},'linear');
        if strcmpi(flag(1:3),'m_m')
            %The FLIPUDs keep M_PATCH from inverting the patch
            h=m_patch(flipud(x1),flipud(y1),str,'edgecolor',str);
        else
            h=patch(x1,y1,str,'edgecolor',str);
        end
    end
end

if ~ishold
    hold off
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function stop = optiplotfval(iter,fval,~)
%OPTIPLOTVAL  Plot objective vs iterations
%
%   Default plotting callback for scalar objectives
%
%   Based on MATLAB's optimplotfval

persistent liter

%Ensure always create a new plot
if(isempty(liter))
    liter = Inf;
end

%Never stop from this function
stop = false;

%Only plot scalars
if(length(fval) > 1)
    fval = norm(fval);
end

%Check for new solver run
if(iter <= liter)
    h = plot(iter,fval,'b*');
    title(sprintf('Current Function Value: %g',fval));
    xlabel('Iteration /  Function Evaluation');
    set(h,'Tag','optiplotf1');
    ylabel('Function Value')
else
    h = findobj(get(gca,'Children'),'Tag','optiplotf1');
    newX = [get(h,'Xdata') iter];
    newY = [get(h,'Ydata') fval];
    set(h,'Xdata',newX, 'Ydata',newY);
    set(get(gca,'Title'),'String',sprintf('Current Function Value: %g',fval));
end

%Save for future comparison
liter = iter;

%Call drawnow to refresh the screen
drawnow;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef RangeAnnotation < handle
    properties (Access = public)
        startSample;
        endSample;
        label;
    end

    properties (Dependent)
        nSamples;
    end

    methods
        function n = get.nSamples(obj)
            n = obj.endSample - obj.startSample + 1;
        end
    end

    methods (Access = public)
        function obj = RangeAnnotation(startSample,endSample,label)
            if nargin > 0
                obj.startSample = startSample;
                obj.endSample = endSample;
                obj.label = label;
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function p =  mcmcProposalProbability2(pE, adjlist, pWhole, smap1, origseg, smap2, newseg, niter)
% compute probability of proposing newseg from smap1
% 1) Randomly select superpixel s1, then randomly selects different superpixel
% s2 within same segment (if one exists); remove s1,s2 from s
% 2) Then, for randomly ordered i:
%     if si is adjacent to s1, assign si to s1 with probability pE(si, s1)
%     if si is adjacent to s2, assign si to s2 with probability pE(si, s2)
%     remove si from s
% 3) Repeat (2) until s is empty

nsp = numel(smap1);

ind1 = find(smap1==origseg);
ind2 = find(smap2==newseg);
n1 = numel(ind1);
n2 = numel(ind2);
if n1 == n2 % entire original segment selected
    if n1==1
        p = 1/nsp;
        disp('one')
    else
        p = pWhole*n2/nsp;
        disp('whole')
    end
else % subset of segment selected
    disp('subset')
    p = (1-pWhole)*n2/nsp;

    adjmat = zeros(n1, n1);
    for k = 1:size(adjlist, 1)
        if smap1(adjlist(k, 1))==origseg && smap1(adjlist(k, 2))==origseg
            i1 = logical(ind1==adjlist(k, 1));
            i2 = logical(ind1==adjlist(k, 2));
            adjmat(i1, i2) = k;
            adjmat(i2, i1) = k;
        end
    end

    assntrue = 2*ones(n1, 1);
    for k = 1:n1
        if any(ind2==ind1(k))
            assntrue(k) = 1;
        end
        adj{k} = find(adjmat(k, :));
        adjlog{k} = logical(adjmat(k, :)>0);
    end
    notselectedind = find(assntrue==2);

    p = p * (n1-n2) / (n1-1); % prob of selecting valid starting point for s2

    cpos = 0;
    cneg = 0;
    assn = zeros(n1, 1);
    for k = 1:niter
        s1 = ind2(ceil(rand(1)*n2));
        ts1 = find(ind1==s1);

        ts2 = notselectedind(ceil(rand(1)*numel(notselectedind)));

        if assntrue(ts2)~=2
            cneg = cneg+1;
        else
            assn(:)=0;
            assn(ts1)=1;
            assn(ts2)=2;
            rind = randperm(n1);
            rind([ts1 ts2]) = [];
            c = 0;
            while ~isempty(rind)
                for t = 1:numel(rind)
                    r = rind(t);
                    if ~assn(r)
                        assnr = assn(adjlog{r});
                    end
                    if ~assn(r) && sum(assnr)>0
                        assnr1 = logical(assnr==1);
                        assnr2 = logical(assnr==2);
                        isadj1 = any(assnr1);
                        isadj2 = any(assnr2);
                        if isadj1 && ~isadj2
                            adj1 = adj{r}(assnr1);
                            p = 1-prod(1-pE(adjmat(r, adj1)));
                            if rand(1) < p
                                assn(r) = 1;
                            end
                        elseif isadj2 && ~isadj1
                            adj2 = adj{r}(assnr2);
                            p = 1-prod(1-pE(adjmat(r, adj2)));
                            if rand(1) < p
                                assn(r) = 2;
                            end
                        else
                            adj1 = adj{r}(assnr1);
                            adj2 = adj{r}(assnr2);
                            p1 = 1-prod(1-pE(adjmat(r, adj1)));
                            p2 = 1-prod(1-pE(adjmat(r, adj2)));
                            p1 = p1 / (p1 + p2);
                            if rand(1) < p1
                                assn(r) = 1;
                            else
                                assn(r) = 2;
                            end
                        end
                        if (assn(r)==1 && assntrue(r)==2) || (assn(r)==2 && assntrue(r)==1)
                            assn(logical(assn==0)) = 3;
                            break;
                        end
                    end % end: if unassigned and has assigned adjacent sp
                end % end: loop through randomly ordered indices
                rind = rind(logical(assn(rind)==0));
            end % end: loop until all assigned while any(~assn)

            if any(assn~=assntrue)
                cneg = cneg + 1;
            else
                cpos = cpos + 1;
            end
        end % end: if ts2 matches profile
    end % end: loop over iterations

    p = p * max(cpos/(cpos+cneg), 1/niter);

end % compute prob for subset

author: zouchuhang

repo: LayoutNet


text:
function table2 = r8mat_border_cut ( m, n, table )

%*****************************************************************************80
%
%% R8MAT_BORDER_CUT cuts the "border" of an R8MAT.
%
%  Discussion:
%
%    We suppose the input data gives values of a quantity on nodes
%    on a 2D grid, and we wish to create a new table corresponding only
%    to those nodes in the interior of the 2D grid.
%
%      0 0 0 0 0 0
%      0 * * * * 0    * * * *
%      0 * * * * 0 -> * * * *
%      0 * * * * 0    * * * *
%      0 0 0 0 0 0
%
%    The illustration suggests the situation in which a 5 by 6 array
%    is input, and a 3 by 4 array is to be output.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the table data.
%
%    Output, real TABLE2(M-2,N-2), the new table data.
%
  if ( m <= 2 || n <= 2 )
    table2 = [];
    return
  end

  table2(1:m-2,1:n-2) = table(2:m-1,2:n-1);

  return
end


text:
%ISCONTOURCONVEX  Tests a contour convexity
%
%     status = cv.isContourConvex(contour)
%
% ## Input
% * __contour__ Input vector of 2D points, stored in numeric array
%   (Nx2/Nx1x2/1xNx2) or cell array of 2-element vectors (`{[x,y], ...}`).
%
% ## Output
% * __status__ Output logical value.
%
% The function tests whether the input contour is convex or not. The contour
% must be simple, that is, without self-intersections. Otherwise, the function
% output is undefined.
%
% See also: cv.convexityDefects
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
report_this_filefun(mfilename('fullpath'));

ttcat = newt2;
mati = maepi(1,3);
%[p,sdp] = mypval(3,mati);
% [p,sdp,c,sdc,dk,sdk,aa,bb]=mypval2(3, mati);

tmin1 = 0.05;

save_aspar2;
do = [ ' ! '  hodi '/aspar/myaspar' ];
eval(do)

load aspar3.out
re = aspar3;

bv = re(2,2);
p =  re(1,2)
c =  re(4,2);
A =  re(3,2);

%[bv magco stan av me mer me2,  pr] =  bvalca3(newt2,1,1);
%[me, bv, si, av] = bmemag(newt2) ;

%A = log10(bv/av)
la = 0;
m0 = maepi(1,6);

m = min(newt2.Magnitude);
dt = 1;

t0 = ( max(newt2.Date) - mati)*365;

c

la = [];ti = [];
for t = c:dt:t0
    la = [la (10^(A + bv*(m0-m)) * (t + c)^(-p))*dt ];
    ti = [ti mati+t/365];
end

figure_w_normalized_uicontrolunits(cum)

l = newt2.Date >= mati + c/365;
tmpn = newt2(l,:);
nu = (1:length(tmpn(:,3))+1); nu(length(tmpn(:,3))+1) = length(tmpn(:,3));

try delete(plc); catch ME, error_handler(ME,@do_nothing);end
try delete(plc2); catch ME, error_handler(ME,@do_nothing);end

hold on
plc = plot([tmpn(:,3) ; teb],nu,'k');

plc2 = plot(ti,cumsum(la),'r');

author: CelsoReyes


text:
Afun2=@(x) (x+1)*(x<1)+(1+1/x)*(x>=1);
fplot(Afun2,[-3,3])

author: Eurus-Holmes


text:
function updatedInfo = matRad_daoVec2ApertureInfo(apertureInfo,apertureInfoVect)
% matRad function to translate vector representation into struct
% The vector representation of the aperture shape and weight are translated
% into an aperture info struct. At the same time, the updated bixel weight
% vector w is computed and a vector listing the correspondence between leaf
% tips and bixel indices for gradient calculation
%
% call
%   updatedInfo = matRad_daoVec2ApertureInfo(apertureInfo,apertureInfoVect)
%
% input
%   apertureInfo:     aperture shape info struct
%   apertureInfoVect: aperture weights and shapes parameterized as vector
%
% output
%   updatedInfo: updated aperture shape info struct according to apertureInfoVect
%
% References
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% function to update the apertureInfo struct after the each iteraton of the
% optimization

w = zeros(apertureInfo.totalNumOfBixels,1);

% initializing variables
updatedInfo = apertureInfo;

updatedInfo.apertureVector = apertureInfoVect;

shapeInd = 1;

indVect = NaN*ones(apertureInfo.totalNumOfShapes + apertureInfo.totalNumOfLeafPairs,1);

% helper function to cope with numerical instabilities through rounding
round2 = @(a,b) round(a*10^b)/10^b;

%% update the shapeMaps
% here the new colimator positions are used to create new shapeMaps that
% now include decimal values instead of binary

% loop over all beams
for i = 1:numel(updatedInfo.beam)

    %posOfRightCornerPixel = apertureInfo.beam(i).posOfCornerBixel(1) + (size(apertureInfo.beam(i).bixelIndMap,2)-1)*apertureInfo.bixelWidth;

    % pre compute left and right bixel edges
    edges_l = updatedInfo.beam(i).posOfCornerBixel(1)...
                 + ([1:size(apertureInfo.beam(i).bixelIndMap,2)]-1-1/2)*updatedInfo.bixelWidth;
    edges_r = updatedInfo.beam(i).posOfCornerBixel(1)...
                 + ([1:size(apertureInfo.beam(i).bixelIndMap,2)]-1+1/2)*updatedInfo.bixelWidth;

    % loop over all shapes
    for j = 1:updatedInfo.beam(i).numOfShapes

        % update the shape weight
        updatedInfo.beam(i).shape(j).weight = apertureInfoVect(shapeInd);

        % get dimensions of 2d matrices that store shape/bixel information
        n = apertureInfo.beam(i).numOfActiveLeafPairs;

        % extract left and right leaf positions from shape vector
        vectorIx     = updatedInfo.beam(i).shape(j).vectorOffset + ([1:n]-1);
        leftLeafPos  = apertureInfoVect(vectorIx);
        rightLeafPos = apertureInfoVect(vectorIx+apertureInfo.totalNumOfLeafPairs);

        % update information in shape structure
        updatedInfo.beam(i).shape(j).leftLeafPos  = leftLeafPos;
        updatedInfo.beam(i).shape(j).rightLeafPos = rightLeafPos;

        % rounding for numerical stability
        leftLeafPos  = round2(leftLeafPos,6);
        rightLeafPos = round2(rightLeafPos,6);

        %
        xPosIndLeftLeaf  = round((leftLeafPos - apertureInfo.beam(i).posOfCornerBixel(1))/apertureInfo.bixelWidth + 1);
        xPosIndRightLeaf = round((rightLeafPos - apertureInfo.beam(i).posOfCornerBixel(1))/apertureInfo.bixelWidth + 1);

        % check limits because of rounding off issues at maximum, i.e.,
        % enforce round(X.5) -> X
        xPosIndLeftLeaf(leftLeafPos == apertureInfo.beam(i).lim_r) = ...
            .5 + (leftLeafPos(leftLeafPos == apertureInfo.beam(i).lim_r) ...
            - apertureInfo.beam(i).posOfCornerBixel(1))/apertureInfo.bixelWidth;
        xPosIndRightLeaf(rightLeafPos == apertureInfo.beam(i).lim_r) = ...
            .5 + (rightLeafPos(rightLeafPos == apertureInfo.beam(i).lim_r) ...
            - apertureInfo.beam(i).posOfCornerBixel(1))/apertureInfo.bixelWidth;

        % find the bixel index that the leaves currently touch
        bixelIndLeftLeaf  = apertureInfo.beam(i).bixelIndMap((xPosIndLeftLeaf-1)*n+[1:n]');
        bixelIndRightLeaf = apertureInfo.beam(i).bixelIndMap((xPosIndRightLeaf-1)*n+[1:n]');

        if any(isnan(bixelIndLeftLeaf)) || any(isnan(bixelIndRightLeaf))
            error('cannot map leaf position to bixel index');
        end

        % store information in index vector for gradient calculation
        indVect(apertureInfo.beam(i).shape(j).vectorOffset+[1:n]-1) = bixelIndLeftLeaf;
        indVect(apertureInfo.beam(i).shape(j).vectorOffset+[1:n]-1+apertureInfo.totalNumOfLeafPairs) = bixelIndRightLeaf;

        % calculate opening fraction for every bixel in shape to construct
        % bixel weight vector

        coveredByLeftLeaf  = bsxfun(@minus,leftLeafPos,edges_l)  / updatedInfo.bixelWidth;
        coveredByRightLeaf = bsxfun(@minus,edges_r,rightLeafPos) / updatedInfo.bixelWidth;

        tempMap = 1 - (coveredByLeftLeaf  + abs(coveredByLeftLeaf))  / 2 ...
                    - (coveredByRightLeaf + abs(coveredByRightLeaf)) / 2;

        % find open bixels
        tempMapIx = tempMap > 0;

        currBixelIx = apertureInfo.beam(i).bixelIndMap(tempMapIx);
        w(currBixelIx) = w(currBixelIx) + tempMap(tempMapIx)*updatedInfo.beam(i).shape(j).weight;

        % save the tempMap (we need to apply a positivity operator !)
        updatedInfo.beam(i).shape(j).shapeMap = (tempMap  + abs(tempMap))  / 2;

        % increment shape index
        shapeInd = shapeInd +1;
    end

end

updatedInfo.bixelWeights = w;
updatedInfo.bixelIndices = indVect;

end

repo: matRad


text:
% fig2i HungarianCV matrix plot

i_fish = 6;
[cIX1,gIX1] = LoadCluster_Direct(i_fish,4,5);
[cIX2,gIX2] = LoadCluster_Direct(i_fish,4,6);

% [score,im1] = HungarianCV(cIX1,cIX2,gIX1,gIX2,isPlotFig);
[score,im1] = HungarianCV(cIX2,cIX1,gIX2,gIX1);

%%
figure;
% subplot(1,2,1)
% imagesc(-im1)
% colormap(bluewhitered)
% axis equal; axis tight;axis xy
%
% subplot(1,2,2)
imagesc(-log(im1))
% colormap(bluewhitered)
axis equal; axis tight;%axis xy
colormap('gray')
ylabel('clusters (1st half)')
xlabel('clusters (2nd half)')
title('Cross-Val: # of cells overlap')
text(size(im1,1)/4,size(im1,2)/20,'score = 0.69');

author: xiuyechen

repo: FishExplorer


text:
function OLPS_cli(dataset)
% OLPS_cli: run OLPS toolbox in CLI (available in both Matlab and Octave)
% this program demos different strategies of on-line portfolio selection
% in the mode of command-line interface (CLI)
%
% OLPS_cli(dataset)
%
% dataset: one can choose any of the following data sets
%  - 'djia':    (default) DJIA (US) perioid 01/14/2001 - 01/14/2003
%  - 'msci':    MSCI (global)       peroid 04/01/2006 - 03/31/2010
%  - 'nyse-n':  NYSE (US)           peroid 01/01/1985 - 06/30/2010
%  - 'nyse-o':  NYSE (US)           peroid 07/03/1962 - 12/31/1984
%  - 'sp500':   S&P500(US)          period 01/02/1998 - 01/31/2003
%  - 'tse':     TSE (CA)            peroid 01/04/1994 - 12/31/1998

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

cd Strategy;

opts.quiet_mode = 1; opts.display_interval = 500;
opts.log_mode = 1; opts.mat_mode = 1;
opts.analyze_mode = 1; opts.progress = 0;
opts.his = 0;

if (nargin<1)
    dataset = 'djia';
end

%----------Benchmarks--------------
manager('ubah', dataset, {0}, opts);
manager('best', dataset, {0}, opts);
manager('ucrp', dataset, {0}, opts);
manager('bcrp', dataset, {0}, opts);

%---------Follow the Winner-------------------
manager('up', dataset, {0}, opts);
manager('eg',  dataset, {0.05, 0}, opts);
manager('ons', dataset, {0, 1, 1/8, 0}, opts);

% manager('sp_start', dataset, {0.25, 0}, opts);
% manager('grw_start', dataset, {0.00005, 0}, opts);
% manager('m0_start', dataset, {0.5, 0}, opts);

%----------Follow the Loser------------------------
manager('anticor', dataset, {30, 0}, opts);
manager('anticor_anticor', dataset, {30, 0}, opts);
manager('pamr', dataset, {0.5, 0}, opts);
manager('pamr_1', dataset, {0.5, 500, 0}, opts);
manager('pamr_2', dataset, {0.5, 500, 0}, opts);
manager('cwmr_var', dataset, {2, 0.5, 0}, opts);
manager('cwmr_stdev', dataset, {2, 0.5, 0}, opts);
manager('olmar1', dataset, {10, 5, 0}, opts);
manager('olmar2', dataset, {10, 0.5, 0}, opts);

%--------Pattern Matching based approach---------
manager('bk', dataset, {5, 10, 1, 0}, opts);
manager('bnn', dataset, {5, 10, 0}, opts);
manager('corn', dataset, {5, 0.1, 0}, opts);
manager('cornu', dataset, {5, 1, 0.1, 0}, opts);
manager('cornk', dataset, {5, 10, 0.1, 0, 1}, opts);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function p = geodesic_create_surface_point(type,id,x,y,z)

p.type = type;
p.id = id;

if nargin == 3
    p.x = x(1);
    p.y = x(2);
    p.z = x(3);
else
    p.x = x;
    p.y = y;
    p.z = z;
end

author: alecjacobson

repo: gptoolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/mmv_phase_transition_snr_40db_s_16.mat';
N = 256;
S = 16;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.SNR = 40;
% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
%% Copyright (C) 2014-2016, 2018-2019, 2022 Colin B. Macdonald
%% Copyright (C) 2014-2015 Andrés Prieto
%% Copyright (C) 2020 Rafael Laboissière
%% Copyright (C) 2020 Jing-Chen Peng
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypemethod  @@sym {@var{sol} =} dsolve (@var{ode})
%% @deftypemethodx @@sym {@var{sol} =} dsolve (@var{ode}, @var{IC})
%% @deftypemethodx @@sym {@var{sol} =} dsolve (@var{ODEs}, @var{IC1}, @var{IC2}, @dots{})
%% @deftypemethodx @@sym {@var{sol} =} dsolve (@var{ODEs}, @var{ICs})
%% @deftypemethodx @@sym {[@var{sol}, @var{classify}] =} dsolve (@dots{})
%% Solve ordinary differential equations (ODEs) symbolically.
%%
%% Basic example:
%% @example
%% @group
%% syms y(x)
%% DE = diff(y, x) - 4*y == 0
%%   @result{} DE = (sym)
%%                 d
%%       -4⋅y(x) + ──(y(x)) = 0
%%                 dx
%% @end group
%%
%% @group
%% sol = dsolve (DE)
%%   @result{} sol = (sym)
%%           4⋅x
%%       C₁⋅ℯ
%% @end group
%% @end example
%%
%% You can specify initial conditions:
%% @example
%% @group
%% sol = dsolve (DE, y(0) == 1)
%%   @result{} sol = (sym)
%%        4⋅x
%%       ℯ
%% @end group
%% @end example
%%
%% In some cases, SymPy can return a classification of the
%% differential equation:
%% @example
%% @group
%% DE = diff(y) == y^2
%%   @result{} DE = (sym)
%%       d           2
%%       ──(y(x)) = y (x)
%%       dx
%%
%% [sol, classify] = dsolve (DE, y(0) == 1)
%%   @result{} sol = (sym)
%%        -1
%%       ─────
%%       x - 1
%%   @result{} classify = ... separable ...
%% @end group
%% @end example
%%
%% Many types of ODEs can be solved, including initial-value
%% problems and boundary-value problem:
%% @example
%% @group
%% DE = diff(y, 2) == -9*y
%%   @result{} DE = (sym)
%%          2
%%         d
%%        ───(y(x)) = -9⋅y(x)
%%          2
%%        dx
%%
%% dsolve (DE, y(0) == 1, diff(y)(0) == 12)
%%   @result{} (sym) 4⋅sin(3⋅x) + cos(3⋅x)
%%
%% dsolve (DE, y(0) == 1, y(sym(pi)/2) == 2)
%%   @result{} (sym) -2⋅sin(3⋅x) + cos(3⋅x)
%% @end group
%% @end example
%%
%% Some systems can be solved, including initial-value problems
%% involving linear systems of first order ODEs with constant
%% coefficients:
%% @example
%% @group
%% syms x(t) y(t)
%% ode_sys = [diff(x(t),t) == 2*y(t);  diff(y(t),t) == 2*x(t)]
%%   @result{} ode_sys = (sym 2×1 matrix)
%%       ⎡d                ⎤
%%       ⎢──(x(t)) = 2⋅y(t)⎥
%%       ⎢dt               ⎥
%%       ⎢                 ⎥
%%       ⎢d                ⎥
%%       ⎢──(y(t)) = 2⋅x(t)⎥
%%       ⎣dt               ⎦
%% @end group
%%
%% @group
%% soln = dsolve (ode_sys)
%%   @result{} soln = scalar structure containing ...
%%        x = ...
%%        y = ...
%%
%% @c doctest: +SKIP_IF(pycall_sympy__ ('return Version(spver) <= Version("1.5.1")'))
%% soln.x
%%   @result{} ans =
%%       (sym)
%%               -2⋅t       2⋅t
%%         - C₁⋅ℯ     + C₂⋅ℯ
%%
%% @c doctest: +SKIP_IF(pycall_sympy__ ('return Version(spver) <= Version("1.5.1")'))
%% soln.y
%%   @result{} ans =
%%       (sym)
%%             -2⋅t       2⋅t
%%         C₁⋅ℯ     + C₂⋅ℯ
%% @end group
%% @end example
%%
%% Note: The Symbolic Math Toolbox used to support strings like 'Dy + y = 0'; we
%% are unlikely to support this so you will need to assemble a symbolic
%% equation instead.
%%
%% @seealso{@@sym/diff, @@sym/int, @@sym/solve}
%% @end deftypemethod

function [soln,classify] = dsolve(ode,varargin)

  % Usually we cast to sym in the _cmd call, but want to be
  % careful here b/c of symfuns
  if (~ iscell (ode) && ~ all (isa (ode, 'sym')))
    error('Inputs must be sym or symfun')
  end

  % FIXME: might be nice to expose SymPy's "sp.ode.classify_sysode" and
  %        "sp.ode.classify_ode" with their own commands
  if (isscalar(ode) && nargout==2)
    cmd = { 'from sympy.solvers import classify_ode'
            'return classify_ode(_ins[0]),' };
    classify = pycall_sympy__ (cmd, ode);
  elseif(~isscalar(ode) && nargout==2)
    warning('Classification of systems of ODEs is currently not supported')
    classify='';
  end

  cmd = { 'ode=_ins[0]; ics=_ins[1:]'
          'if len(ics) == 1:'
          '    ics = ics[0]'
          'try:'
          '    ics = {ic.lhs: ic.rhs for ic in ics}'
          'except TypeError:'  % not iterable
          '    ics = {ics.lhs: ics.rhs}'
          'sol = sp.dsolve(ode, ics=ics)'
          'def convert_helper(sympy_obj):'
          '    if isinstance(sympy_obj, Eq) and sympy_obj.lhs.is_Function:'
                   % y(t) = rhs to str "y", rhs expression
          '        return str(sympy_obj.lhs.func), sympy_obj.rhs'
          '    return None, None'
          % if just single equality with simple lhs, we return the rhs
          'if isinstance(sol, Eq):'
          '    if sol.lhs.is_Function:'
          '        return sol.rhs'
          % If the solution set is iterable (system or multiple solutions),
          % we will try to convert to structure of {"x": expr, "y": expr2, ...}
          'try:'
          '    return_data = dict()'
          '    for solution_part in sol:'
          '        key, rhs = convert_helper(solution_part)'
          '        if key is None:'
          '            raise ValueError("not of right form for extraction")'
          '        if key in return_data.keys():'
          '            raise KeyError(f"repeated key {key}")'
          '        return_data[key] = rhs'
          '    return return_data'
          'except (TypeError, ValueError, KeyError):'
          '    pass'
          % if nothing else worked, give back whatever form we have
          'return sol' };
  soln = pycall_sympy__ (cmd, ode, varargin{:});
end

%!error <sym> dsolve (1, sym('x'))

%!test
%! syms y(x)
%! de = diff(y, 2) - 4*y == 0;
%! f = dsolve(de);
%! syms C1 C2
%! g1 = C1*exp(-2*x) + C2*exp(2*x);
%! g2 = C2*exp(-2*x) + C1*exp(2*x);
%! assert (isequal (f, g1) || isequal (f, g2))

%!test
%! % Not enough initial conditions
%! syms y(x) C1
%! de = diff(y, 2) + 4*y == 0;
%! g = 3*cos(2*x) + C1*sin(2*x);
%! try
%!   f = dsolve(de, y(0) == 3);
%!   waserr = false;
%! catch
%!   waserr = true;
%!   expectederr = regexp (lasterr (), 'Perhaps.*under-specified');
%!   f = 42;
%! end
%! assert ((waserr && expectederr) || isequal (f, g))

%!test
%! % Solution in implicit form
%! syms y(x) C1
%! de = (2*x*y(x) - exp(-2*y(x)))*diff(y(x), x) + y(x) == 0;
%! sol = dsolve (de);
%! eqn = x*exp(2*y(x)) - log(y(x)) == C1;
%! % could differ by signs
%! sol = lhs (sol) - rhs (sol);
%! eqn = lhs (eqn) - rhs (eqn);
%! sol2 = subs (sol, C1, -C1);
%! assert (isequal (sol, eqn) || isequal (sol2, eqn))

%%!xtest
%%! % system with solution in implicit form
%%! % TODO: not implemented upstream?
%%! syms y(x) z(x) C1
%%! de1 = (2*x*y(x) - exp(-2*y(x)))*diff(y(x), x) + y(x) == 0;
%%! de2 = diff(z, x) == 0;
%%! sol = dsolve ([de1; de2]);

%!test
%! % Compute solution and classification
%! syms y(x) C1
%! de = (2*x*y(x) - exp(-2*y(x)))*diff(y(x), x) + y(x) == 0;
%! [sol, classy] = dsolve (de);
%! assert (any (strcmp (classy, '1st_exact')))

%!test
%! % initial conditions (first order ode)
%! syms y(x)
%! de = diff(y, x) + 4*y == 0;
%! f = dsolve(de, y(0) == 3);
%! g = 3*exp(-4*x);
%! assert (isequal (f, g))

%!test
%! % initial conditions (second order ode)
%! syms y(x)
%! de = diff(y, 2) + 4*y == 0;
%! f = dsolve(de, y(0) == 3, subs(diff(y,x),x,0)==0);
%! g = 3*cos(2*x);
%! assert (isequal (f, g))

%!test
%! % Dirichlet boundary conditions (second order ode)
%! syms y(x)
%! de = diff(y, 2) + 4*y == 0;
%! f = dsolve(de, y(0) == 2, y(1) == 0);
%! g = -2*sin(2*x)/tan(sym('2'))+2*cos(2*x);
%! assert (isequal (simplify (f - g), 0))

%!test
%! % Neumann boundary conditions (second order ode)
%! syms y(x)
%! de = diff(y, 2) + 4*y == 0;
%! f = dsolve(de, subs(diff(y,x),x,0)==1, subs(diff(y,x),x,1)==0);
%! g = sin(2*x)/2+cos(2*x)/(2*tan(sym('2')));
%! assert (isequal (simplify (f - g), 0))

%!test
%! % Dirichlet-Neumann boundary conditions (second order ode)
%! syms y(x)
%! de = diff(y, 2) + 4*y == 0;
%! f = dsolve(de, y(0) == 3, subs(diff(y,x),x,1)==0);
%! g = 3*sin(2*x)*tan(sym('2'))+3*cos(2*x);
%! assert (isequal (simplify (f - g), 0))

%!test
%! % System of ODEs gives struct, Issue #1003.
%! syms x(t) y(t)
%! ode1 = diff(x(t),t) == 2*y(t);
%! ode2 = diff(y(t),t) == 2*x(t);
%! soln = dsolve([ode1, ode2]);
%! assert (isstruct (soln))
%! assert (numfields (soln) == 2)
%! assert (isequal (sort (fieldnames (soln)), {'x'; 'y'}))

%!test
%! % System of ODEs
%! syms x(t) y(t) C1 C2
%! ode1 = diff(x(t),t) == 2*y(t);
%! ode2 = diff(y(t),t) == 2*x(t);
%! soln = dsolve([ode1, ode2]);
%! soln = [soln.x, soln.y];
%! g1 = [C1*exp(-2*t) + C2*exp(2*t), -C1*exp(-2*t) + C2*exp(2*t)];
%! g2 = [C1*exp(2*t) + C2*exp(-2*t), C1*exp(2*t) - C2*exp(-2*t)];
%! g3 = [-C1*exp(-2*t) + C2*exp(2*t), C1*exp(-2*t) + C2*exp(2*t)];
%! g4 = [C1*exp(2*t) - C2*exp(-2*t), C1*exp(2*t) + C2*exp(-2*t)];
%! % old SymPy <= 1.5.1 had some extra twos
%! g5 = [2*C1*exp(-2*t) + 2*C2*exp(2*t), -2*C1*exp(-2*t) + 2*C2*exp(2*t)];
%! g6 = [2*C1*exp(2*t) + 2*C2*exp(-2*t), 2*C1*exp(2*t) - 2*C2*exp(-2*t)];
%! assert (isequal (soln, g1) || isequal (soln, g2) || ...
%!         isequal (soln, g3) || isequal (soln, g4) || ...
%!         isequal (soln, g5) || isequal (soln, g6))

%!test
%! % System of ODEs (initial-value problem)
%! syms x(t) y(t)
%! ode_1=diff(x(t),t) == 2*y(t);
%! ode_2=diff(y(t),t) == 2*x(t);
%! sol_ivp=dsolve([ode_1,ode_2],x(0)==1,y(0)==0);
%! g_ivp=[exp(-2*t)/2+exp(2*t)/2,-exp(-2*t)/2+exp(2*t)/2];
%! assert (isequal ([sol_ivp.x, sol_ivp.y], g_ivp))

%!test
%! syms y(x)
%! de = diff(y, 2) + 4*y == 0;
%! f = dsolve(de, y(0) == 0, y(sym(pi)/4) == 1);
%! g = sin(2*x);
%! assert (isequal (f, g))

%!test
%! % Nonlinear example
%! syms y(x) C1
%! e = diff(y, x) == y^2;
%! g = -1 / (C1 + x);
%! soln = dsolve(e);
%! assert (isequal (soln, g))

%!test
%! % Nonlinear example with initial condition
%! syms y(x)
%! e = diff(y, x) == y^2;
%! g = -1 / (x - 1);
%! soln = dsolve(e, y(0) == 1);
%! assert (isequal (soln, g))

%!test
%! % forcing, Issue #183, broken in older sympy
%! if (pycall_sympy__ ('return Version(spver) >= Version("1.7.1")'))
%! syms x(t) y(t)
%! ode1 = diff(x) == x + sin(t) + 2;
%! ode2 = diff(y) == y - t - 3;
%! soln = dsolve([ode1 ode2], x(0) == 1, y(0) == 2);
%! X = soln.x;
%! Y = soln.y;
%! assert (isequal (diff(X) - (X + sin(t) + 2), 0))
%! assert (isequal (diff(Y) - (Y - t - 3), 0))
%! end

%!test
%! syms f(x) a b
%! de = diff(f, x) == 4*f;
%! s = dsolve(de, f(a) == b);
%! assert (isequal (subs(s, x, a), b))

%!test
%! % array of ICs
%! syms x(t) y(t)
%! ode_1 = diff (x(t), t) == 2*y(t);
%! ode_2 = diff (y(t), t) == 2*x(t);
%! sol = dsolve([ode_1, ode_2], [x(0)==1 y(0)==0]);
%! g = [exp(-2*t)/2+exp(2*t)/2, -exp(-2*t)/2+exp(2*t)/2];
%! assert (isequal ([sol.x, sol.y], g))

%!test
%! % cell-array of ICs or ODEs, but not both
%! % Note: to support both we'd need a wrapper outside of @sym
%! syms x(t) y(t)
%! ode_1 = diff (x(t), t) == 2*y(t);
%! ode_2 = diff (y(t), t) == 2*x(t);
%! sol = dsolve([ode_1, ode_2], {x(0)==1 y(0)==0});
%! g = [exp(-2*t)/2+exp(2*t)/2, -exp(-2*t)/2+exp(2*t)/2];
%! assert (isequal ([sol.x, sol.y], g))
%! sol = dsolve({ode_1, ode_2}, [x(0)==1 y(0)==0]);
%! g = [exp(-2*t)/2+exp(2*t)/2, -exp(-2*t)/2+exp(2*t)/2];
%! assert (isequal ([sol.x, sol.y], g))

%!test
%! % array of ICs, Issue #1040.
%! if (pycall_sympy__ ('return Version(spver) >= Version("1.7.1")'))
%! syms x(t) y(t) z(t)
%! syms x_0 y_0 z_0
%! diffEqns = [diff(x, t) == -x + 1, diff(y, t) == -y, diff(z, t) == -z];
%! initCond = [x(0) == x_0, y(0) == y_0, z(0) == z_0];
%! soln = dsolve (diffEqns, initCond);
%! soln = [soln.x, soln.y, soln.z];
%! exact_soln = [(x_0 - 1)*exp(-t) + 1  y_0*exp(-t)  z_0*exp(-t)];
%! assert (isequal (soln, exact_soln))
%! end

repo: octsympy


text:
function setupDataOpts(vocYear, testName, datasetDir)

global DATAopts;

% Setup VOC data
devkitroot = [datasetDir, 'VOCdevkit', '/'];
DATAopts.year = vocYear;
DATAopts.dataset = sprintf('VOC%d', DATAopts.year);
DATAopts.datadir        = [devkitroot, DATAopts.dataset, '/'];
DATAopts.resdir         = [devkitroot, 'results', '/', DATAopts.dataset '/'];
DATAopts.localdir       = [devkitroot, 'local', '/', DATAopts.dataset, '/'];
DATAopts.gStructPath    = [DATAopts.resdir, 'GStructs', '/'];
DATAopts.imgsetpath     = [DATAopts.datadir, 'ImageSets', '/', 'Main', '/', '%s.txt'];
DATAopts.imgpath        = [DATAopts.datadir, 'JPEGImages', '/', '%s.jpg'];
DATAopts.clsimgsetpath  = [DATAopts.datadir, 'ImageSets', '/', 'Main', '/', '%s_%s.txt'];
DATAopts.annopath       = [DATAopts.datadir, 'Annotations', '/', '%s.xml'];
DATAopts.annocachepath	= [DATAopts.localdir, '%s_anno.mat'];
DATAopts.classes={...
    'aeroplane'
    'bicycle'
    'bird'
    'boat'
    'bottle'
    'bus'
    'car'
    'cat'
    'chair'
    'cow'
    'diningtable'
    'dog'
    'horse'
    'motorbike'
    'person'
    'pottedplant'
    'sheep'
    'sofa'
    'train'
    'tvmonitor'};
DATAopts.nclasses = length(DATAopts.classes);
DATAopts.testset = testName;
DATAopts.minoverlap = 0.5;

author: nightrome

repo: matconvnet-calvin


text:
function a = sin(a)
%SIN          Gradient sine sin(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  ax = cos(full(a.x(:)));
  a.x = sin(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function sparse_grid_open_test05 ( dim_num, level_max )

%*****************************************************************************80
%
%% TEST05 tests LEVELS_OPEN_INDEX to make a Gauss Patterson grid.
%
%  Discussion:
%
%    This routine gets the sparse grid indices and determines the
%    corresponding sparse grid abscissas for a Gauss Patterson scheme.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05:\n' );
  fprintf ( 1, '  Make a sparse Gauss Patterson grid.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVELS_OPEN_INDEX returns all grid indexes\n' );
  fprintf ( 1, '  whose level value satisfies\n' );
  fprintf ( 1, '    0 <= LEVEL <= LEVEL_MAX.\n' );
  fprintf ( 1, '  Here, LEVEL is the sum of the levels of the 1D rules,\n' );
  fprintf ( 1, '  and the order of the rule is 2^(LEVEL+1) - 1.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now we demonstrate how to convert grid indices\n' );
  fprintf ( 1, '  into physical grid points.  In this case, we\n' );
  fprintf ( 1, '  want points on [-1,+1]^DIM_NUM.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_MAX = %d\n', level_max );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );

  point_num = sparse_grid_ofn_size ( dim_num, level_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of unique points in the grid = %d\n', point_num );
%
%  Compute the orders and points.
%
  grid_index = levels_open_index ( dim_num, level_max, point_num );
%
%  Now we're done.  Print the merged grid data.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Grid index:\n' );
  fprintf ( 1, '\n' );
  for j = 1 : point_num
    fprintf ( 1, '  %4d', j );
    for dim = 1 : dim_num
      fprintf ( 1, '%6d', grid_index(dim,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Convert index information to physical information.
%
  grid_point = zeros ( dim_num, point_num );

  order_max = 2^( level_max + 1 ) - 1;

  for j = 1 : point_num
    for dim = 1 : dim_num
      grid_point(dim,j) = gp_abscissa ( order_max, grid_index(dim,j) );
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Grid points:\n' );
  fprintf ( 1, '\n' );
  for j = 1 : point_num
    fprintf ( 1, '  %4d', j );
    for dim = 1 : dim_num
      fprintf ( 1, '%10f', grid_point(dim,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestDrawAxis
    %TestDrawAxis

    methods (Static)
        function test_1
            img = zeros([50 50 3], 'uint8');
            camMat = eye(3);
            distCoef = zeros(1,5);
            rvec = zeros(1,3);
            tvec = zeros(1,3);
            len = 1;
            out = cv.drawAxis(img, camMat, distCoef, rvec, tvec, len);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.drawAxis();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end


text:
function gb=gabor_fn(sigma,theta,lambda,psi,gamma)

sigma_x = sigma;
sigma_y = sigma/gamma;

% Bounding box
nstds = 3;
xmax = max(abs(nstds*sigma_x*cos(theta)),abs(nstds*sigma_y*sin(theta)));
xmax = ceil(max(1,xmax));
ymax = max(abs(nstds*sigma_x*sin(theta)),abs(nstds*sigma_y*cos(theta)));
ymax = ceil(max(1,ymax));
xmin = -xmax; ymin = -ymax;
[x,y] = meshgrid(xmin:xmax,ymin:ymax);

% Rotation
x_theta=x*cos(theta)+y*sin(theta);
y_theta=-x*sin(theta)+y*cos(theta);

gb= 1/(2*pi*sigma_x *sigma_y) * exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+psi);

author: Sable

repo: mcbench-benchmarks


text:
function fd1d_advection_diffusion_steady ( )

%*****************************************************************************80
%
%% FD1D_ADVECTION_DIFFUSION_STEADY solves steady advection diffusion equation.
%
%  Discussion:
%
%    The steady advection diffusion equation has the form:
%
%      v ux - k * uxx = 0
%
%    where V (the advection velocity) and K (the diffusivity) are positive
%    constants, posed in the region
%
%      a = 0 < x < 1 = b
%
%    with boundary conditions
%
%      u(0) = 0, u(1) = 1.
%
%    The discrete solution is unreliable when dx > 2 * k / v / ( b - a ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Hall, Thomas Porsching,
%    Numerical Analysis of Partial Differential Equations,
%    Prentice-Hall, 1990,
%    ISBN: 013626557X,
%    LC: QA374.H29.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_ADVECTION_DIFFUSION_STEADY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the 1D steady advection diffusion equation:\n' );
  fprintf ( 1, '    v du/dx - k d2u/dx2 = 0\n' );
  fprintf ( 1, '  with constant, positive velocity V and diffusivity K,\n' );
  fprintf ( 1, '  over the interval:\n' );
  fprintf ( 1, '    0.0 <= x <= 1.0\n' );
  fprintf ( 1, '  with boundary conditions:\n' );
  fprintf ( 1, '    u(0) = 0, u(1) = 1.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use finite differences:\n' );
  fprintf ( 1, '    d u/dx  = (u(x+dx)-u(x-dx))/2/dx\n' );
  fprintf ( 1, '    d2u/dx2 = (u(x+dx)-2u(x)+u(x-dx))/dx^2\n' );
%
%  Physical constants.
%
  v = 1.0;
  k = 0.05;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Diffusivity K = %g\n', k );
  fprintf ( 1, '  Constant velocity V = %g\n', v );
%
%  Spatial discretization.
%
  nx = 41;
  a = 0.0;
  b = 1.0;
  dx = ( b - a ) / ( nx - 1 );
  x = ( linspace ( a, b, nx ) )';

  fprintf ( 1, '  Number of nodes NX = %d\n', nx );
  fprintf ( 1, '  DX = %g\n', dx );
  fprintf ( 1, '  Maximum safe DX is %g\n', 2 * k / v / ( b - a ) );
%
%  Set up linear system corresponding to boundary conditions and
%  advection-diffusion equation.
%
  A = sparse ( nx, nx );
  f = zeros ( nx, 1 );

  A(1,1) = 1.0;
  f(1) = 0.0;

  for i = 2 : nx - 1
    A(i,i-1) = - v / dx / 2.0 -       k / dx / dx;
    A(i,i)   =                + 2.0 * k / dx / dx;
    A(i,i+1) = + v / dx / 2.0 -       k / dx / dx;
    f(i) = 0.0;
  end

  A(nx,nx) = 1.0;
  f(nx) = 1.0;
%
%  Solve linear system for U.
%
  u = A \ f;
%
%  The exact solution to the differential equation is known.
%
  r = v * ( b - a ) / k;
  w = ( 1.0 - exp ( r * x ) ) / ( 1.0 - exp ( r ) );
%
%  Plot discrete solution U (dots) and exact solution W (line).
%
  hold on
  plot ( x, w, 'b-', 'LineWidth', 2 );
  plot ( x, u, 'r.', 'LineWidth', 2, 'Markersize', 20 );

  grid on
  xlabel ( '<--X-->', 'Fontsize', 16 );
  ylabel ( '<--U(X)-->', 'Fontsize', 16 );
  title ( 'Solution of Steady Advection-Diffusion Equation', 'Fontsize', 24);
  hold off

  filename = 'fd1d_advection_diffusion_steady.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot as "%s".\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_ADVECTION_DIFFUSION_STEADY\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function Y = fft2_mid0(X)
% fft2_mid0 -- 2d fft with argument [-pi,pi] (instead of [0,2pi])
% Usage:
%   Y = fft_mid0(X)
% Inputs:
%   X	Array(n)
% Outputs:
%   Y   Array(n)
% Description:
%  Performs 1d fft with grid (-n/2):(n/2-1) on both time
%  and frequency side.
%
% y(k) = sum_{-n/2 <= t1,t2 < n/2} exp(-i 2pi/n (k1 t1 + k2 t2)) x(t1,t2),
%                                         (-n/2) <= k1,k2 < n/2
%

Y = fftshift(fft2(fftshift(X)));

author: thomaskuestner


text:
function [ssimval, ssimmap] = ssim(varargin)
%SSIM Structural Similarity Index for measuring image quality
%   SSIMVAL = SSIM(A, REF) calculates the Structural Similarity Index
%   (SSIM) value for image A, with the image REF as the reference. A and
%   REF can be 2D grayscale or 3D volume images, and must be of the same
%   size and class.
%
%   [SSIMVAL, SSIMMAP] = SSIM(A, REF) also returns the local SSIM value for
%   each pixel in SSIMMAP. SSIMMAP has the same size as A.
%
%   [SSIMVAL, SSIMMAP] = SSIM(A, REF, NAME1, VAL1,...) calculates the SSIM
%   value using name-value pairs to control aspects of the computation.
%   Parameter names can be abbreviated.
%
%   Parameters include:
%
%   'Radius'                 - Specifies the standard deviation of
%                              isotropic Gaussian function used for
%                              weighting the neighborhood pixels around a
%                              pixel for estimating local statistics. This
%                              weighting is used to avoid blocking
%                              artifacts in estimating local statistics.
%                              The default value is 1.5.
%
%   'DynamicRange'           - Positive scalar, L, that specifies the
%                              dynamic range of the input image. By
%                              default, L is chosen based on the class of
%                              the input image A, as L =
%                              diff(getrangefromclass(A)). Note that when
%                              class of A is single or double, L = 1 by
%                              default.
%
%   'RegularizationConstants'- Three-element vector, [C1 C2 C3], of
%                              non-negative real numbers that specifies the
%                              regularization constants for the luminance,
%                              contrast, and structural terms (see [1]),
%                              respectively. The regularization constants
%                              are used to avoid instability for image
%                              regions where the local mean or standard
%                              deviation is close to zero. Therefore, small
%                              non-zero values should be used for these
%                              constants. By default, C1 = (0.01*L).^2, C2
%                              = (0.03*L).^2, and C3 = C2/2, where L is the
%                              specified 'DynamicRange' value. If a value
%                              of 'DynamicRange' is not specified, the
%                              default value is used (see name-value pair
%                              'DynamicRange').
%
%   'Exponents'               - Three-element vector [alpha beta gamma],
%                               of non-negative real numbers that specifies
%                               the exponents for the luminance, contrast,
%                               and structural terms (see [1]),
%                               respectively. By default, all the three
%                               exponents are 1, i.e. the vector is [1 1
%                               1].
%
%   Notes
%   -----
%   1. A and REF can be arrays of upto three dimensions. All 3D arrays
%      are considered 3D volumetric images. RGB images will also be
%      processed as 3D volumetric images.
%
%   2. Input image A and reference image REF are converted to
%      floating-point type for internal computation.
%
%   3. For signed-integer images (int16), an offset is applied to bring the
%      gray values in the non-negative range before computing the SSIM
%      index.
%
%   Example
%   ---------
%   This example shows how to compute SSIM value for a blurred image given
%   the original reference image.
%
%   ref = imread('pout.tif');
%   H = fspecial('Gaussian',[11 11],1.5);
%   A = imfilter(ref,H,'replicate');
%
%   subplot(1,2,1); imshow(ref); title('Reference Image');
%   subplot(1,2,2); imshow(A);   title('Blurred Image');
%
%   [ssimval, ssimmap] = ssim(A,ref);
%
%   fprintf('The SSIM value is %0.4f.\n',ssimval);
%
%   figure, imshow(ssimmap,[]);
%   title(sprintf('SSIM Index Map - Mean SSIM Value is %0.4f',ssimval));
%
%   Class Support
%   -------------
%   Input arrays A and REF must be one of the following classes: uint8,
%   int16, uint16, single, or double. Both A and REF must be of the same
%   class. They must be nonsparse. SSIMVAL is a scalar and SSIMMAP is an
%   array of the same size as A. Both SSIMVAL and SSIMMAP are of class
%   double, unless A and REF are of class single in which case SSIMVAL and
%   SSIMMAP are of class single.
%
%   References:
%   -----------
%   [1] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
%       Quality Assessment: From Error Visibility to Structural
%       Similarity," IEEE Transactions on Image Processing, Volume 13,
%       Issue 4, pp. 600- 612, 2004.
%
%   See also IMMSE, MEAN, MEDIAN, PSNR, SUM, VAR.

%   Copyright 2013-2014 The MathWorks, Inc.

narginchk(2,10);

[A, ref, C, exponents, radius] = parse_inputs(varargin{:});

if isempty(A)
    ssimval = zeros(0, 'like', A);
    ssimmap = A;
    return;
end

if isa(A,'int16') % int16 is the only allowed signed-integer type for A and ref.
    % Add offset for signed-integer types to bring values in the
    % non-negative range.
    A = double(A) - double(intmin('int16'));
    ref = double(ref) - double(intmin('int16'));
elseif isinteger(A)
    A = double(A);
    ref = double(ref);
end

% Gaussian weighting function
gaussFilt = getGaussianWeightingFilter(radius,ndims(A));

% Weighted-mean and weighted-variance computations
mux2 = imfilter(A, gaussFilt,'conv','replicate');
muy2 = imfilter(ref, gaussFilt,'conv','replicate');
muxy = mux2.*muy2;
mux2 = mux2.^2;
muy2 = muy2.^2;

sigmax2 = imfilter(A.^2,gaussFilt,'conv','replicate') - mux2;
sigmay2 = imfilter(ref.^2,gaussFilt,'conv','replicate') - muy2;
sigmaxy = imfilter(A.*ref,gaussFilt,'conv','replicate') - muxy;

% Compute SSIM index
if (C(3) == C(2)/2) && isequal(exponents(:),ones(3,1))
    % Special case: Equation 13 from [1]
    num = (2*muxy + C(1)).*(2*sigmaxy + C(2));
    den = (mux2 + muy2 + C(1)).*(sigmax2 + sigmay2 + C(2));
    if (C(1) > 0) && (C(2) > 0)
        ssimmap = num./den;
    else
        % Need to guard against divide-by-zero if either C(1) or C(2) is 0.
        isDenNonZero = (den ~= 0);
        ssimmap = ones(size(A));
        ssimmap(isDenNonZero) = num(isDenNonZero)./den(isDenNonZero);
    end

else
    % General case: Equation 12 from [1]
    % Luminance term
    if (exponents(1) > 0)
        num = 2*muxy + C(1);
        den = mux2 + muy2 + C(1);
        ssimmap = guardedDivideAndExponent(num,den,C(1),exponents(1));
    else
        ssimmap = ones(size(A), 'like', A);
    end

    % Contrast term
    sigmaxsigmay = [];
    if (exponents(2) > 0)
        sigmaxsigmay = sqrt(sigmax2.*sigmay2);
        num = 2*sigmaxsigmay + C(2);
        den = sigmax2 + sigmay2 + C(2);
        ssimmap = ssimmap.*guardedDivideAndExponent(num,den,C(2),exponents(2));
    end

    % Structure term
    if (exponents(3) > 0)
        num = sigmaxy + C(3);
        if isempty(sigmaxsigmay)
            sigmaxsigmay = sqrt(sigmax2.*sigmay2);
        end
        den = sigmaxsigmay + C(3);
        ssimmap = ssimmap.*guardedDivideAndExponent(num,den,C(3),exponents(3));
    end

end

ssimval = mean(ssimmap(:));

end

% -------------------------------------------------------------------------
function component = guardedDivideAndExponent(num, den, C, exponent)

if C > 0
    component = num./den;
else
    component = ones(size(num),'like',num);
    isDenNonZero = (den ~= 0);
    component(isDenNonZero) = num(isDenNonZero)./den(isDenNonZero);
end

if (exponent ~= 1)
    component = component.^exponent;
end

end

function gaussFilt = getGaussianWeightingFilter(radius,N)
% Get 2D or 3D Gaussian weighting filter

filtRadius = ceil(radius*3); % 3 Standard deviations include >99% of the area.
filtSize = 2*filtRadius + 1;

if (N < 3)
    % 2D Gaussian mask can be used for filtering even one-dimensional
    % signals using imfilter.
    gaussFilt = fspecial('gaussian',[filtSize filtSize],radius);
else
    % 3D Gaussian mask
     [x,y,z] = ndgrid(-filtRadius:filtRadius,-filtRadius:filtRadius, ...
                    -filtRadius:filtRadius);
     arg = -(x.*x + y.*y + z.*z)/(2*radius*radius);
     gaussFilt = exp(arg);
     gaussFilt(gaussFilt<eps*max(gaussFilt(:))) = 0;
     sumFilt = sum(gaussFilt(:));
     if (sumFilt ~= 0)
         gaussFilt  = gaussFilt/sumFilt;
     end
end

end

function [A, ref, C, exponents, radius] = parse_inputs(varargin)

validImageTypes = {'uint8','uint16','int16','single','double'};

A = varargin{1};
validateattributes(A,validImageTypes,{'nonsparse','real'},mfilename,'A',1);

ref = varargin{2};
validateattributes(ref,validImageTypes,{'nonsparse','real'},mfilename,'REF',2);

if ~isa(A,class(ref))
    error(message('images:validate:differentClassMatrices','A','REF'));
end

if ~isequal(size(A),size(ref))
    error(message('images:validate:unequalSizeMatrices','A','REF'));
end

if (ndims(A) > 3)
    error(message('images:validate:tooManyDimensions','A and REF',3));
end

% Default values for parameters
dynmRange = diff(getrangefromclass(A));
C = [];
exponents = [1 1 1];
radius = 1.5;

args_names = {'dynamicrange', 'regularizationconstants','exponents',...
              'radius'};

for i = 3:2:nargin
    arg = varargin{i};
    if ischar(arg)
        idx = find(strncmpi(arg, args_names, numel(arg)));
        if isempty(idx)
            error(message('images:validate:unknownInputString', arg))

        elseif numel(idx) > 1
            error(message('images:validate:ambiguousInputString', arg))

        elseif numel(idx) == 1
            if (i+1 > nargin)
                error(message('images:validate:missingParameterValue'));
            end
            if idx == 1
                dynmRange = varargin{i+1};
                validateattributes(dynmRange,{'numeric'},{'positive', ...
                    'finite', 'real', 'nonempty','scalar'}, mfilename, ...
                    'DynamicRange',i);
                dynmRange = double(dynmRange);

            elseif idx == 2
                C = varargin{i+1};
                validateattributes(C,{'numeric'},{'nonnegative','finite', ...
                    'real','nonempty','vector', 'numel', 3}, mfilename, ...
                    'RegularizationConstants',i);
                C = double(C);

            elseif idx == 3
                exponents = varargin{i+1};
                validateattributes(exponents,{'numeric'},{'nonnegative', ...
                    'finite', 'real', 'nonempty','vector', 'numel', 3}, ...
                    mfilename,'Exponents',i);
                exponents = double(exponents);

            elseif idx == 4
                radius = varargin{i+1};
                validateattributes(radius,{'numeric'},{'positive','finite', ...
                    'real', 'nonempty','scalar'}, mfilename,'Radius',i);
                radius = double(radius);
            end
        end
    else
        error(message('images:validate:mustBeString'));
    end
end

% If 'RegularizationConstants' is not specified, choose default C.
if isempty(C)
    C = [(0.01*dynmRange).^2 (0.03*dynmRange).^2 ((0.03*dynmRange).^2)/2];
end

end

author: rwenqi

repo: GFN-dehazing

save_path: github-repos/MATLAB/rwenqi-GFN-dehazing


text:
% function [nll, grad] = InstanceNegLogLikelihood(X, y, theta, modelParams)
% returns the negative log-likelihood and its gradient, given a CRF with parameters theta,
% on data (X, y).
%
% Inputs:
% X            Data.                           (numCharacters x numImageFeatures matrix)
%              X(:,1) is all ones, i.e., it encodes the intercept/bias term.
% y            Data labels.                    (numCharacters x 1 vector)
% theta        CRF weights/parameters.         (numParams x 1 vector)
%              These are shared among the various singleton / pairwise features.
% modelParams  Struct with three fields:
%   .numHiddenStates     in our case, set to 26 (26 possible characters)
%   .numObservedStates   in our case, set to 2  (each pixel is either on or off)
%   .lambda              the regularization parameter lambda
%
% Outputs:
% nll          Negative log-likelihood of the data.    (scalar)
% grad         Gradient of nll with respect to theta   (numParams x 1 vector)
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function [nll, grad] = InstanceNegLogLikelihood(X, y, theta, modelParams)

    % featureSet is a struct with two fields:
    %    .numParams - the number of parameters in the CRF (this is not numImageFeatures
    %                 nor numFeatures, because of parameter sharing)
    %    .features  - an array comprising the features in the CRF.
    %
    % Each feature is a binary indicator variable, represented by a struct
    % with three fields:
    %    .var          - a vector containing the variables in the scope of this feature
    %    .assignment   - the assignment that this indicator variable corresponds to
    %    .paramIdx     - the index in theta that this feature corresponds to
    %
    % For example, if we have:
    %
    %   feature = struct('var', [2 3], 'assignment', [5 6], 'paramIdx', 8);
    %
    % then feature is an indicator function over X_2 and X_3, which takes on a value of 1
    % if X_2 = 5 and X_3 = 6 (which would be 'e' and 'f'), and 0 otherwise.
    % Its contribution to the log-likelihood would be theta(8) if it's 1, and 0 otherwise.
    %
    % If you're interested in the implementation details of CRFs,
    % feel free to read through GenerateAllFeatures.m and the functions it calls!
    % For the purposes of this assignment, though, you don't
    % have to understand how this code works. (It's complicated.)

    featureSet = GenerateAllFeatures(X, modelParams);

    % Use the featureSet to calculate nll and grad.
    % This is the main part of the assignment, and it is very tricky - be careful!
    % You might want to code up your own numerical gradient checker to make sure
    % your answers are correct.
    %
    % Hint: you can use CliqueTreeCalibrate to calculate logZ effectively.
    %       We have halfway-modified CliqueTreeCalibrate; complete our implementation
    %       if you want to use it to compute logZ.

    nll = 0;
    grad = zeros(size(theta));
    %%%
    % Your code here:
    n = length(y);
    K = modelParams.numHiddenStates;
    Factors = repmat(struct('var',[],'card',[],'val',[]),1,2*n-1);
    for i = 1:n-1
	    Factors(i).var = [i];
	    Factors(i).card = [K];
	    Factors(i).val = zeros(1,prod(Factors(i).card));
	    Factors(n+i).var = [i,i+1];
	    Factors(n+i).card = [K,K];
	    Factors(n+i).val = zeros(1,prod(Factors(i+n).card));
    end
    Factors(n).var = [n];
    Factors(n).card = [K];
    Factors(n).val = zeros(1,prod(Factors(n).card));

    featureCounts = zeros(size(theta));
    modelFeatureCounts = zeros(size(theta));

    for i = 1:length(featureSet.features)
	  feature = featureSet.features(i);
	  if length(feature.var) == 1
		  Factors(feature.var(1)).val(feature.assignment) += theta(feature.paramIdx);
	  else
		  % assuming feature.var is always [i,i+1].. if this fails accomodate for [i+1,i] as well
	          indx = AssignmentToIndex(feature.assignment,[K,K]);
		  fact = n + feature.var(1);
		  Factors(fact).val(indx) = theta(feature.paramIdx);
	  end

	  if feature.assignment == y([feature.var])
		  featureCounts(feature.paramIdx) = 1;
	  end
    end

    weightedFeatureCounts = featureCounts.*theta;

    logFactors = Factors;

    for i = 1:length(Factors)
	    Factors(i).val = exp(Factors(i).val);
    end

    P = CreateCliqueTree(Factors);
    [P, logZ] = CliqueTreeCalibrate(P, 0);
    lambda = modelParams.lambda;
    regularizationCost = sum(theta.*theta)*lambda/2;
    regularizationGradient = theta*lambda;

    % change the code how we do inference over all the probablilities to in case of words with letters not equal to 3.. do it using code instead of this lame crap..

    p1 = FactorMarginalization(P.cliqueList(1),2);
    p2 = FactorMarginalization(P.cliqueList(1),1);
    p3 = FactorMarginalization(P.cliqueList(2),2);
    p12 = P.cliqueList(1);
    p23 = P.cliqueList(2);
    p1.val/=sum(p1.val);
    p2.val/=sum(p2.val);
    p3.val/=sum(p3.val);
    p23.val/=sum(p23.val);
    p12.val/=sum(p12.val);
    ft = pppp(featureSet.features,p1,p2,p3,p12,p23);
    for i = 1:length(ft)
	    modelFeatureCounts(ft(i).paramIdx) += ft(i).p;
    end

    nll = logZ - sum(weightedFeatureCounts) + regularizationCost;
    grad = modelFeatureCounts - featureCounts + regularizationGradient;
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function M = mult(disc, f)
%MULT    Multiplication operator for the TRIGSPEC class.
%   M = MULT(A, F) returns the multiplication operator that represents
%   u(x) -> F(x)u(x), in the Fourier basis.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Obtaining some useful information:
n = disc.dimension;
dom = disc.domain;

% Convert to a TRIGTECH-based CHEBFUN. Use a non-adaptive construction:
% we know that the TRIGTECH-based CHEBFUN will have a length smaller
% than the CHEBTECH-based CHEBFUN.
f = chebfun(f, dom, 'periodic', length(f));
f = simplify(f);

% Call TRIGSPEC/MULTMAT.
M = trigspec.multmat(n, f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function dist_signed = triangle_point_dist_signed_2d ( t, p )

%*****************************************************************************80
%
%% TRIANGLE_POINT_DIST_SIGNED_2D: signed distance ( triangle, point ) in 2D.
%
%  Discussion:
%
%    If the signed distance is:
%    0, the point is on the boundary of the triangle;
%    negative, the point is in the triangle;
%    positive, the point is outside the triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%    These should be given in counter clockwise order.
%
%    Input, real P(2,1), the point which is to be checked.
%
%    Output, real DIST_SIGNED, the signed distance from the
%    point to the triangle.
%
  dim_num = 2;
%
%  Compute the signed line distances to the point.
%
  dis12 = line_exp_point_dist_signed_2d ( t(1:2,1), t(1:2,2), p );

  dis23 = line_exp_point_dist_signed_2d ( t(1:2,2), t(1:2,3), p );

  dis31 = line_exp_point_dist_signed_2d ( t(1:2,3), t(1:2,1), p );
%
%  If the point is inside the triangle, all the line distances are negative.
%  The largest (negative) line distance has the smallest magnitude,
%  and is the signed triangle distance.
%
  if ( dis12 <= 0.0 & dis23 <= 0.0 & dis31 <= 0.0 )

    dist_signed = max ( dis12, max ( dis23, dis31 ) );
%
%  If the point is outside the triangle, then we have to compute
%  the (positive) line segment distances and take the minimum.
%
  else

    dis12 = segment_point_dist_2d ( t(1:2,1), t(1:2,2), p );
    dis23 = segment_point_dist_2d ( t(1:2,2), t(1:2,3), p );
    dis31 = segment_point_dist_2d ( t(1:2,3), t(1:2,1), p );

    dist_signed = min ( dis12, min ( dis23, dis31 ) );

  end

  return
end