text:
function [W, funcVal] = Least_NCFGLassoF2(X, Y, lambda, beta, opts)
% Non-Convex Fused Group Lasso with Least Squares Loss (Formulation 2).
% By Jiayu Zhou (jiayu.zhou@asu.edu) Dec. 2011

% Objective (Non-Convex):
% argmin_W { \sum_{i=1}^t (0.5 * norm (Y{i} - X{i}' * W(:, i))^2)
%    + lambda * \sum_{i=1}^d \sqrt{ \|R * W(i, :)\|_1 + beta  * \|W(i, :)\|_1}
%          }
%
% t - task number
% d - dimension
% n - sample size
%
% lambda: joint feature selection
% beta:   Lasso sparsity.
% R encodes fused structure relationship [1 -1 0 ...; 0 1 -1 ...; ...]
%    R=zeros(t,t-1);R(1:(t+1):end)=1;R(2:(t+1):end)=-1;

%%% input
% X: {d * n} * t - input matrix
% Y: {n * 1} * t - output matrix
% lambda: 1        - joint feature selection
% beta:   1        - Lasso sparsity.

if nargin <5
    opts = [];
end

if isfield(opts, 'max_iter')
    max_iter = opts.max_iter;
else
    max_iter = 20;
end

if isfield(opts, 'tol_funcVal')
    tol_funcVal = opts.tol_funcVal;
else
    tol_funcVal = 10^-6;
end

task_num  = length (X);
dimension = size(X{1}, 1);

% Relation
R=zeros(task_num,task_num-1);
R(1:(task_num+1):end)=1;
R(2:(task_num+1):end)=-1;
R = R';

if isfield(opts, 'starting_point')
    %fprintf('Least_NCFGLassoF2: Init model used\n');
    W0 = opts.starting_point;
else
    %W0 = rand(dimension, task_num);
    [W0] = Least_FGLasso(X, Y, lambda*beta, lambda);
end

epsilon = 10e-12; % should be slight larger than machine error.

funcVal = funcVal_eval (W0);

inner_iteration_opts.max_iter = 100;
inner_iteration_opts.tol_funcVal = 10^-7;

iter = 0;
while iter < max_iter
    % reweigting vector.

    rho1 = zeros(dimension, 1); % weighted Lasso
    rho2 = zeros(dimension, 1); % weighted Fused Lasso
    for t = 1: dimension
        rho2(t) = lambda / (2 * sqrt(norm(W0(t, :),1) + norm(R * W0(t, :)', 1) + epsilon));
        rho1(t) = rho2(t) * beta;
    end

    [W, sub_funcVal] = Least_Weight2FGLasso(X, Y, rho1, rho2, inner_iteration_opts);

    funcVal = cat(1, funcVal, funcVal_eval (W));

    % test stop condition.
    if length(funcVal)> 1 &&...
            (abs(funcVal(end-1) - funcVal(end))/ funcVal(end))<=tol_funcVal
        break;
    end

    iter = iter + 1;
    W0 = W;

end

    function [funcVal] = funcVal_eval (W)
        funcVal = 0;
        for i = 1: task_num
            funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
        end
        % non-smooth part
        for i = 1 : size(W, 1) % dimension
            w = W(i, :);
            funcVal = funcVal ...
                + lambda * sqrt(norm(R * w', 1) + beta * norm(w, 1));
        end
    end

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
classdef HuberLoss < dagnn.Loss

  properties
    sigma = 1.
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnhuberloss(inputs{1}, inputs{2}, ...
                       'instanceWeights', inputs{3}, ...
                       obj.opts{:}) ;
      n = obj.numAveraged ;
      m = n + size(inputs{1}, 4) ;
      obj.average = (n * obj.average + gather(outputs{1})) / m ;
      obj.numAveraged = m ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} = vl_nnhuberloss(inputs{1}, inputs{2}, derOutputs{1}, ...
                                'instanceWeights', inputs{3}, ...
                                obj.opts{:}) ;
      derInputs{2} = [] ;
      derInputs{3} = [] ;
      derParams = {} ;
    end

    function obj = HuberLoss(varargin)
      obj.load(varargin) ;
      obj.sigma = obj.sigma ;
    end
  end
end

repo: MFT


text:
function [Pw] = RFT_Pval(u,k,c,fwhm,L,type,dof)
% computes Pval according to Friston et al. (1996), Eq. 1-4.
% [Pw] = RFT_Pval(u,k,c,fwhm,L)
% This function ca be called to derive the p-value at different levels of
% inference, i.e.:
% - P_peak = RFT_Pval(u,0,1,fwhm,L)
% - P_cluster = RFT_Pval(U,k,1,fwhm,L), where U was the set-inducing
% threshold, and k is the observed spatial extent of the cluster.
% - P_set = RFT_Pval(U,K,c,fwhm,L), where U and K were the set-inducing
% thresholds, and c was the number of observed upcrossing clusters
% IN:
%   - u: RF value
%   - k: cluster's spatial extent
%   - c: number of upcrossing clusters
%   - fwhm: estimated FWHM
%   - L: search volume
%   - type: type of RF. Can be set to 'norm' (normal, default), 't'
%   (Student) or 'F' (Fisher).
%   - dof: degrees of freedom (only relevant for 't' or 'F' fields).
% OUT:
%   - Pw: the ensuing p-value

EC = RFT_expectedTopo(u,L,fwhm,1,type,dof);
switch type
    case 'norm'
        P0 = 1-VBA_spm_Ncdf(u,0,1);
    case 't'
        P0 = 1-VBA_spm_Tcdf(u,dof);
    case 'F'
        P0 = 1-VBA_spm_Fcdf(u,dof(1),dof(2));
end
beta = (gamma(3/2).*EC./(L.*P0)).^2;
Pnk = exp(-beta.*k.^2);
if k == 0, Pnk = 1; end; % solves the numerical issue when P0=0
Pw = 1;
for i=0:c-1
    Pw = Pw - myPoissonPMF(i,EC.*Pnk);
end

function p = myPoissonPMF(x,Ex)
p = (Ex.^x).*exp(-Ex)./factorial(x);

% function p = myPoissonCDF(x,Ex)
% p = 1 - gammainc(Ex,x+1);

author: MBB-team

repo: VBA-toolbox


text:
classdef prtRegress < prtAction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prtRegress is an abstract base class for all regression objects.

    properties
        plotOptions = prtRegress.initializePlotOptions(); % Plotting Options
    end

    properties (SetAccess = protected)
        isSupervised = true; % True
        isCrossValidateValid = true; % True
	end

    methods (Hidden = true)
        function featureNameModificationFunction = getFeatureNameModificationFunction(obj)
			featureNameModificationFunction = prtUtilFeatureNameModificationFunctionHandleCreator('%s Output_{#index#}', obj.nameAbbreviation);
        end
	end

    methods
        function obj = prtRegress()
            % As an action subclass we must set the properties to reflect
            % our dataset requirements
            obj.classTrain = 'prtDataSetRegress';
            obj.classRun = 'prtDataSetStandard';
            obj.classRunRetained = true;
        end

        function varargout = plot(Obj)
            % PLOT  Plot the prtRegress object
            %
            % OBJ.plot() plots a trained prtRegress object. The plot
            % displays the original data points, the regressed data points,
            % and a line or curve interpolating between the regressed data
            % points.

            assert(Obj.isTrained,'Regressor must be trained before it can be plotted.');

            [OutputDataSet, linGrid] = runRegressorOnGrid(Obj);

            colors = Obj.plotOptions.colorsFunction(Obj.dataSetSummary.nTargetDimensions);
            lineWidth = Obj.plotOptions.lineWidth;

            switch Obj.dataSetSummary.nFeatures
                case 1
                    HandleStructure.regressorPlotHandle = plot(linGrid,OutputDataSet.getObservations,'color',colors(1,:),'lineWidth',lineWidth);
                case 2
                    sz = ones(1,Obj.dataSetSummary.nFeatures)*Obj.plotOptions.nSamplesPerDim(Obj.dataSetSummary.nFeatures);
                    HandleStructure.regressorPlotHandle = surf(...
                        reshape(linGrid(:,1),sz),...
                        reshape(linGrid(:,2),sz),...
                        reshape(OutputDataSet.getObservations,sz));
                    HandleStructure.regressorPlotHandle.EdgeColor = colors(1,:);
                    HandleStructure.regressorPlotHandle.FaceColor = 'none';
                    HandleStructure.regressorPlotHandle.LineWidth = lineWidth;
                otherwise
                    error('nFeatures in the training dataset must be 1 or 2');
            end

            holdState = get(gca,'nextPlot');
            if ~isempty(Obj.dataSet)
                hold on
                HandleStructure.dataSetPlotHandle = plot(Obj.dataSet);
            end
            set(gca,'nextPlot',holdState);

            axis tight;
            title(Obj.name)

            varargout = {};
            if nargout > 0
                varargout = {HandleStructure};
            end
        end

        function [OutputDataSet, linGrid, gridSize] = runRegressorOnGrid(Obj, upperBounds, lowerBounds)

            if nargin < 3 || isempty(lowerBounds)
                lowerBounds = Obj.dataSetSummary.lowerBounds;
            end

            if nargin < 2 || isempty(upperBounds)
                upperBounds = Obj.dataSetSummary.upperBounds;
            end

            [linGrid, gridSize] = prtPlotUtilGenerateGrid(upperBounds, lowerBounds, Obj.plotOptions.nSamplesPerDim);

            OutputDataSet = run(Obj,prtDataSetRegress(linGrid));
        end
    end
    methods (Static, Hidden = true)
        function plotOptions = initializePlotOptions()
            if ~isdeployed
                plotOptions = prtOptionsGet('prtOptionsRegressPlot');
            else
                plotOptions = prtOptions.prtOptionsRegressPlot;
            end
        end
    end
end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% Copyright (C) 2010 Quan Wang <wangq10@rpi.edu>
% Signal Analysis and Machine Perception Laboratory
% Department of Electrical, Computer, and Systems Engineering
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA

%% A demo showing how to use this package for binary classification
%    f: (n0+n1)*2 feature matrix, each row being a data point
%    l0: (n0+n1)*1 ground truth of binary label vector, each element being 0 or 1
%    l: (n0+n1)*1 resulting binary label vector, each element being 0 or 1
clear;clc;close all;

n0=100;
n1=100;

l0=[zeros(n0,1);ones(n1,1)];
f=[l0+randn(n0+n1,1),l0+randn(n0+n1,1)];

[w,t,fp]=fisher_training(f,l0);
[l,precision,recall,accuracy,F1]=fisher_testing(f,w,t,l0);

%% visualization
figure;
plot(f(1:n0,1),f(1:n0,2),'bo','MarkerSize',10);
hold on;
plot(f(n0+1:end,1),f(n0+1:end,2),'rs','MarkerSize',10);
grid on;

xx=-2:0.1:3;
yy=-w(1)/w(2)*xx+t/w(2);
plot(xx,yy,'-.k','LineWidth',2);
title('data points and classification border');
legend('label 0','label 1','class border');

author: Sable

repo: mcbench-benchmarks


text:
function r = struct2rowsexport(s,varargin)
% r = struct2rowsexport(s,varargin)
% Optional PVs:
%   - nansAsStrings: If true, numeric scalar NaNs are replaced with the
%   string 'NaN'. Defaults to false.

nansAsStrings = myparse(varargin,...
  'nansAsStrings',false);

assert(isstruct(s)&&isvector(s),'Expected d to be a struct vector.');
s = s(:);

r = struct2cell(s);
r = r';

% replace values in r for export
for c = 1:numel(r)
  val = r{c};
  if isscalar(val)&&isnan(val)&&nansAsStrings
    r{c} = 'NaN';
  elseif isnumeric(val)&&isscalar(val) || ...
      islogical(val)&&isscalar(val) || ...
      ischar(val)&&isrow(val) || ...
      isequal(val,[])
    % none
  else
    r{c} = '<array value>';
  end
end

hdr = fieldnames(s);
r = [hdr(:)';r];

author: kristinbranson

repo: JAABA


text:
function [f,err] = fmincon_funn(cbData,nRow,x,njdiff,dXjbase,reserved,inParam)

persistent params hash xevaled oldh F
if nargin == 7
    params = inParam;
    hash = randn(1,length(params.linearindicies));
    oldh = randn(1);
    xevaled = zeros(1,length(params.interfacedata.c));
    F = [];
    return
end

x = x(1:length(params.linearindicies));

if x(:)'*hash(:) ~= oldh
    oldh = x(:)'*hash(:);
    xevaled = zeros(1,length(params.interfacedata.c));
    xevaled(params.linearindicies) = x;

    % Experimental support for arbitrary functions
    if ~isempty(params.interfacedata.evalMap)

        pp = params.monomtable(params.nonlinearindicies,:);
        xevaled(params.nonlinearindicies) = prod(repmat(xevaled,length(params.nonlinearindicies),1).^pp,2);

        for i = 1:length(params.interfacedata.evalMap)
            arguments = {params.interfacedata.evalMap{i}.fcn,xevaled(params.interfacedata.evalMap{i}.variableIndex)};
            arguments = {arguments{:},params.interfacedata.evalMap{i}.arg{2:end-1}};
            xevaled(params.interfacedata.evalVariables(i)) = feval(arguments{:});
        end
    end
    xevaled(params.nonlinearindicies) = prod(repmat(xevaled,length(params.nonlinearindicies),1).^params.monomtable(params.nonlinearindicies,:),2);
    xevaled = xevaled(:);
    F =  -params.F_struc*[1;xevaled];
end

if nRow == -1
    f = params.interfacedata.c'*xevaled+xevaled'*params.interfacedata.Q*xevaled;
else
    f = F(nRow + 1);%-params.F_struc(nRow + 1,:)*[1;xevaled];
end
err = 0;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [graph_prop, graph_prop_glob] = bct_toolbox_undirected_graph_metrics(r, thresh, varargin)
% Calculate some Sporns BCT toolbox functions
%
% graph_prop = bct_toolbox_undirected_graph_metrics(r, [threshold_input])
%
% Inputs:
% r = correlation matrix
% thresh = 0 to 1 value. (.1 is a common value)
%
% optional:
%   'doplots' - show matrices
%   'doweighted' -- this increase compute time by several orders of
%   magnitude. Default = false
%
% Outputs:
% graph_prop = A table of node-level graph metrics
% graph_prop_glob = A table of global graph metrics
%
% For BCT toolbox, see:
% https://sites.google.com/site/bctnet/
%
% For descriptions of metrics, see:
% https://sites.google.com/site/bctnet/measures/list
%
% Notes:
% - Some metrics are calculated on thresholded, weighted (continuous valued) matrix r
%   Sig matrix or link density threshold will be applied.
%   To use without any thresholding, sig mat should be all ones, or link density 100
% - Others are calculated on thresholded, binarized matrix (positive connections only)
%   If you enter a sig matrix, this will be used as the input to graph metric functions
%   If you enter a link density, a binary matrix will be calculated
% - Fisher r to Z is computed, which will have little impact in many cases
%   but may help in some cases
% - Many BCT functions will not use negative values and ignore weights, so require thresholding.
%   P-values/statistical significance is one way to threshold.
%   Another way would be using multiple arbitrary thresholds (e.g., 10% link density)
%   The latter may be practical for individual subjects, unless individual-subject stats are calculated and saved
% - This function prioritizes
% - Degree, etc. may possibly be related to confounds (head motion), and may want to adjust/test for this.
%
% Examples:
% -------------------------------------------------------
% % start with r, a series of correlation matrices, one per subject (k x k x n_subjects)
% OUT = ttest3d(r);                                             % Get some basic stats
% graph_prop = bct_toolbox_undirected_graph_metrics(OUT.r);     % input group mean correlation matrix
%
% rr = corr(table2array(graph_prop));                           % correlate the metrics
% plot_correlation_matrix(rr, 'names', graph_prop.Properties.VariableNames)
%
% % Note: correlations among community vectors are not meaningful.
% % for similarity among communities/modules detected with different methods, see partition_distance.m

% Check path
if isempty(which('threshold_proportional.m'))
    error('The BCT toolbox does not seem to be on your Matlab path. See https://sites.google.com/site/bctnet/');
end

doplots = false;
doweighted = false;
if any(strcmp(varargin, 'doplots')), doplots = true; end
if any(strcmp(varargin, 'doweighted')), doweighted = true; end

% Prep r for BCT undirected
r = double(r);
r = (r' + r) ./ 2;          % enforce symmetry (rounding error possible)
r = r - eye(size(r));       % for BCT and squareform

% Missing regions/constant values will give NaNs, so need to account for them

numnans = any(isnan(r));
if sum(numnans) > .10 * length(numnans)
    warning('More than 10% of vars have NaN correlation values [missing data?]. Replacing NaNs with 0s. Be careful about effects on subsequent metrics');
end

r(isnan(r)) = 0;

% Threshold: Use sig matrix or link density
bu_matrix = weight_conversion(threshold_proportional(r, thresh), 'binarize');

% Make weighted matrix based on Fisher transform
z = rToZ(r);
z(isinf(z)) = max(z(~isinf(z))); % correlations of 1.00 get transformed to Inf, which causes BCT to crash in some functions. Replace w/ max value. Note that this code does not correctly handle case of r = -1.00, as this is very unlikely.
wu_matrix = z;

% view
if doplots
    create_figure('BCT networks', 1, 2)
    imagesc(bu_matrix), colorbar
    subplot(1,2,2)
    imagesc(wu_matrix), colorbar
end

% Node-level properties: Undirected binary and weighted networks
graph_prop = table();

graph_prop.community = modularity_und(bu_matrix, 1);           % C: communities from adjacency matrix
graph_prop.core_w = core_periphery_dir(wu_matrix, 1)';         % Core-periphery

graph_prop.degree = degrees_und(bu_matrix)';                   % Node degree
graph_prop.betweenness_bin = betweenness_bin(double(bu_matrix))';  % note: logical did not work in some cases...
graph_prop.clustercoef_bin = clustering_coef_bu(bu_matrix);        % Clustering coefficient
graph_prop.local_efficiency_bin = efficiency_bin(bu_matrix, 1);    % Local efficiency

% weighted matrix computations
if doweighted
    graph_prop.strength = strengths_und(wu_matrix)';               % weighted strength. tested: same results as using r_mean
    graph_prop.local_efficiency_weighted = efficiency_wei(wu_matrix, 1);    % Local efficiency
    graph_prop.eigenvector_centrality = eigenvector_centrality_und(wu_matrix);    % Eigenvector centrality (similar to PageRank)
    graph_prop.clustercoef_weighted = clustering_coef_wu(wu_matrix);        % Clustering coefficient
    graph_prop.betweenness_weighted = betweenness_wei(wu_matrix);  % note: slow

    graph_prop.community_w = modularity_und(wu_matrix, 1);
    graph_prop.community_w2 = community_louvain(wu_matrix, 1, [], 'negative_asym');
end
    % Others to consider
% assortativity_bin(bu_network, 0);
% kcore_bu
% graph_prop.betweenness_w = betweenness_wei(OUT.r);            % slow
% graph_prop.rich_club_coeff = rich_club_bu(bu_network)';       % output len not same as number of nodes

% Note: for similarity among communities/modules, see partition_distance.m

% Global properties
graph_prop_glob = table();

graph_prop_glob.trans_bu = transitivity_bu(bu_matrix);
graph_prop_glob.glob_efficiency = efficiency_bin(bu_matrix);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function boundary = p11_boundary_nearest ( m, n, point )

%*****************************************************************************80
%
%% P11_BOUNDARY_NEAREST returns a nearest boundary point in problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, real BOUNDARY(M,N), points on the boundary
%    that are nearest to each point.
%
  q1 = [ ...
    0.00, 0.00; ...
    1.00, 0.00; ...
    0.75, 0.25; ...
    0.25, 0.25 ]';
  q2 = [ ...
    0.00, 0.00; ...
    0.25, 0.25; ...
    0.25, 0.75; ...
    0.00, 1.00 ]';
  q3 = [ ...
    0.50, 0.50; ...
    0.50, 0.25; ...
    0.75, 0.25; ...
    1.00, 0.50 ]';
  q4 = [ ...
    0.50, 0.50; ...
    0.25, 0.50; ...
    0.25, 0.25; ...
    0.50, 0.25 ]';
  q5 = [ ...
    0.50, 0.50; ...
    0.50, 1.00; ...
    0.25, 0.75; ...
    0.25, 0.50 ]';
  t1 = [ ...
    1.00, 0.00; ...
    1.00, 0.50; ...
    0.75, 0.25 ]';
  t2 = [ ...
    0.00, 1.00; ...
    0.25, 0.25; ...
    0.50, 1.00 ]';
  t3 = [ ...
    0.50, 0.50; ...
    1.00, 0.50; ...
    1.00, 1.00 ]';
  t4 = [ ...
    0.50, 0.50; ...
    1.00, 1.00; ...
    0.50, 1.00 ]';

  x1 =  0.0;
  x2 =  0.5;
  x3 = +1.0;
  y1 =  0.0;
  y2 =  0.5;
  y3 = +1.0;

  for j = 1 : n

    if ( point(1,j) <= x1 & point(2,j) <= y1 )

      boundary(1,j) = x1;
      boundary(2,j) = y1;

    elseif ( point(1,j) <= x1 & point(2,j) <= y3 )

      boundary(1,j) = x1;
      boundary(2,j) = point(2,j);

    elseif ( point(1,j) <= x1 & y3 <= point(2,j) )

      boundary(1,j) = x1;
      boundary(2,j) = y3;

    elseif ( point(1,j) <= x3 & point(2,j) <= y1 )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y1;

    elseif ( quad_contains_point_2d ( q1, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y1;

    elseif ( quad_contains_point_2d ( q2, point(1:2,j) ) )

      boundary(1,j) = x1;
      boundary(2,j) = point(2,j);

    elseif ( triangle_contains_point_2d ( t1, point(1:2,j) ) )

      boundary(1,j) = x3;
      boundary(2,j) = point(2,j);

    elseif ( quad_contains_point_2d ( q3, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y2;

    elseif ( quad_contains_point_2d ( q4, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = y2;

    elseif ( quad_contains_point_2d ( q5, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = point(2,j);

    elseif ( triangle_contains_point_2d ( t2, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y3;

    elseif ( x1 <= point(1,j) & point(1,j) <= x2 & y3 <= point(2,j) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y3;

    elseif ( triangle_contains_point_2d ( t3, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y2;

    elseif ( triangle_contains_point_2d ( t4, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = point(2,j);

    elseif ( x2 <= point(1,j) & point(1,j) <= x3 & y3 <= point(2,j) )

      boundary(1,j) = x2;
      boundary(2,j) = y3;

    elseif ( x3 <= point(1,j) & point(2,j) <= y1 )

      boundary(1,j) = x3;
      boundary(2,j) = y1;

    elseif ( x3 <= point(1,j) & y1 <= point(2,j) & point(2,j) <= y2 )

      boundary(1,j) = x3;
      boundary(2,j) = point(2,j);

    elseif ( x3 <= point(1,j) & y2 <= point(2,j) & point(2,j) <= y3 )

      boundary(1,j) = x3;
      boundary(2,j) = y2;

    elseif ( x3 <= point(1,j) & y3 <= point(2,j) & point(2,j) <= point(1,j) )

      boundary(1,j) = x3;
      boundary(2,j) = y2;

    elseif ( x3 <= point(1,j) & y3 <= point(2,j) & point(1,j) <= point(2,j) )

      boundary(1,j) = x2;
      boundary(2,j) = y3;

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'P11_BOUNDARY_NEAREST - Fatal error!\n' );
      fprintf ( 1, '  Orphan point = ( %f, %f )\n', point(1:2,j) );
      error ( 'P11_BOUNDARY_NEAREST - Fatal error!' );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function amanatidesWooAlgorithm(origin, direction, grid3D, verbose)
% A fast and simple voxel traversal algorithm through a 3D space partition (grid)
% proposed by J. Amanatides and A. Woo (1987).
%
% Input:
%    origin.
%    direction.
%    grid3D: grid dimensions (nx, ny, nz, minBound, maxBound).
% Author:
%    Jesús P. Mena-Chalco.

    if (verbose)
        figure;
        hold on;
        text(origin(1), origin(2), origin(3), 'origin');
        plot3(origin(1), origin(2), origin(3), 'k.', 'MarkerSize', 15);
        quiver3(origin(1), origin(2), origin(3), direction(1), direction(2), direction(3), 30);

        vmin = grid3D.minBound';
        vmax = grid3D.maxBound';
        BoxVertices = [vmax(1) vmin(2) vmin(3); vmax(1) vmax(2) vmin(3); vmin(1) vmax(2) vmin(3); vmin(1) vmax(2) vmax(3); vmin(1) vmin(2) vmax(3); vmax(1) vmin(2) vmax(3); vmin; vmax ];
        BoxFaces = [1 2 3 7; 1 2 8 6; 1 6 5 7; 7 5 4 3; 2 8 4 3; 8 6 5 4];
        h = patch('Vertices',BoxVertices,'Faces',BoxFaces,'FaceColor','yellow');
        set(h, 'FaceAlpha', 0.1);

        view(60,30);
        axis tight;
        xlabel('x');
        ylabel('y');
        zlabel('z');
        grid on;
    end;

    [flag, tmin] = rayBoxIntersection(origin, direction, grid3D.minBound, grid3D.maxBound);

    if (flag==0)
        disp('\n The ray does not intersect the grid');
    else
        if (tmin<0)
            tmin = 0;
        end;

        start   = origin + tmin*direction;
        boxSize = grid3D.maxBound-grid3D.minBound;

        if (verbose)
            plot3(start(1), start(2), start(3), 'r.', 'MarkerSize', 15);
        end;

        x = floor( ((start(1)-grid3D.minBound(1))/boxSize(1))*grid3D.nx )+1;
        y = floor( ((start(2)-grid3D.minBound(2))/boxSize(2))*grid3D.ny )+1;
        z = floor( ((start(3)-grid3D.minBound(3))/boxSize(3))*grid3D.nz )+1;

        if (x==(grid3D.nx+1));  x=x-1;  end;
        if (y==(grid3D.ny+1));  y=y-1;  end;
        if (z==(grid3D.nz+1));  z=z-1;  end;

        if (direction(1)>=0)
            tVoxelX = (x)/grid3D.nx;
            stepX = 1;
        else
            tVoxelX = (x-1)/grid3D.nx;
            stepX = -1;
        end;

        if (direction(2)>=0)
            tVoxelY = (y)/grid3D.ny;
            stepY = 1;
        else
            tVoxelY = (y-1)/grid3D.ny;
            stepY = -1;
        end;

        if (direction(3)>=0)
            tVoxelZ = (z)/grid3D.nz;
            stepZ = 1;
        else
            tVoxelZ = (z-1)/grid3D.nz;
            stepZ = -1;
        end;

        voxelMaxX  = grid3D.minBound(1) + tVoxelX*boxSize(1);
        voxelMaxY  = grid3D.minBound(2) + tVoxelY*boxSize(2);
        voxelMaxZ  = grid3D.minBound(3) + tVoxelZ*boxSize(3);

        tMaxX      = tmin + (voxelMaxX-start(1))/direction(1);
        tMaxY      = tmin + (voxelMaxY-start(2))/direction(2);
        tMaxZ      = tmin + (voxelMaxZ-start(3))/direction(3);

        voxelSizeX = boxSize(1)/grid3D.nx;
        voxelSizeY = boxSize(2)/grid3D.ny;
        voxelSizeZ = boxSize(3)/grid3D.nz;

        tDeltaX    = voxelSizeX/abs(direction(1));
        tDeltaY    = voxelSizeY/abs(direction(2));
        tDeltaZ    = voxelSizeZ/abs(direction(3));

        while ( (x<=grid3D.nx)&&(x>=1) && (y<=grid3D.ny)&&(y>=1) && (z<=grid3D.nz)&&(z>=1) )

            if (verbose)
                fprintf('\nIntersection: voxel = [%d %d %d]', [x y z]);

                t1 = [(x-1)/grid3D.nx, (y-1)/grid3D.ny, (z-1)/grid3D.nz ]';
                t2 = [  (x)/grid3D.nx,  (y)/grid3D.ny,    (z)/grid3D.nz ]';

                vmin = (grid3D.minBound + t1.*boxSize)';
                vmax = (grid3D.minBound + t2.*boxSize)';

                smallBoxVertices = [vmax(1) vmin(2) vmin(3); vmax(1) vmax(2) vmin(3); vmin(1) vmax(2) vmin(3); vmin(1) vmax(2) vmax(3); vmin(1) vmin(2) vmax(3); vmax(1) vmin(2) vmax(3); vmin; vmax ];
                smallBoxFaces    = [1 2 3 7; 1 2 8 6; 1 6 5 7; 7 5 4 3; 2 8 4 3; 8 6 5 4];

                h = patch('Vertices', smallBoxVertices, 'Faces', smallBoxFaces, 'FaceColor', 'blue', 'EdgeColor', 'white');
                set(h,'FaceAlpha',0.2);
            end;

            % ---------------------------------------------------------- %
            % check if voxel [x,y,z] contains any intersection with the ray
            %
            %   if ( intersection )
            %       break;
            %   end;
            % ---------------------------------------------------------- %

            if (tMaxX < tMaxY)
                if (tMaxX < tMaxZ)
                    x = x + stepX;
                    tMaxX = tMaxX + tDeltaX;
                else
                    z = z + stepZ;
                    tMaxZ = tMaxZ + tDeltaZ;
                end;
            else
                if (tMaxY < tMaxZ)
                    y = y + stepY;
                    tMaxY = tMaxY + tDeltaY;
                else
                    z = z + stepZ;
                    tMaxZ = tMaxZ + tDeltaZ;
                end;
            end;
        end;
     end;
end


text:
function planC = insertDeformS(planC,deformS)

indexS = planC{end};

deformIndex = length(planC{indexS.deform}) + 1;
planC{indexS.deform}  = dissimilarInsert(planC{indexS.deform},deformS,deformIndex);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function y=log10p1(x)
%%LOG10P1 Evaluate the function log10(x+1) without a loss of precision for
%         values of x near 0.
%
%INPUTS: x A vector of matrix of real values at which the function should
%          be evaluated.
%
%OUTPUTS: y The value(s) of the function log10(x+1) at the point(s) in x.
%
%For values of x having magnitude over 1e-2 and for complex values, the
%function is evaluated as written. For smaller values, a Taylor series
%expansion is used to avoid the singularity near the origin.
%
%EXAMPLE:
%Compare
% log10(1e-16+1)
% log10p1(1e-16)
%For the first one, one will get 0, for the second,
%4.342944819032518e-17. Compare this to a value computed using
%higher-precision arithmetic of 4.3429448190325180593640482375402*10^-17
%and it is clear that this function is more accurate than explicitly
%evaluating the expression.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Allocate space
y=zeros(size(x));

%For large enough values, use the standard function.
sel=abs(x)>1e-2|~isreal(x);
x1=x(sel);
y(sel)=log10(x1+1);

%For smaller magnitude values, use an ninth-order Taylor series expansion.
%The relative error at 10^-2 is on the order of 1e-20, which is less than
%the precision of a double floating point number at that value. The series
%is given in Horner form.
log10Val=log(10);
x1=x(~sel);
y(~sel)=x1.*(x1.*(x1.*(x1.*(x1.*(x1.*(x1.*(x1.*(-(1/(8*log10Val))+x/(9*log10Val))+1/(7*log10Val))-1/(6*log10Val))+1/(5*log10Val))-1/(4*log10Val))+1/(3*log10Val))-1/(2*log10Val))+1/log10Val);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% housekeeping
close all
clear
clc
%% add the necessary paths
rise_startup()

%% we read the model
m=dsge('fs2000');
%% we solve, passing a function that solves the steady state analytically

mwssf_a=solve(m,'steady_state_file','fs2000_steadystate');

%% we solve passing a function that gives initial values of the steady state.
% this is the counterpart to dynare's initval

mwssf_b=solve(m,'steady_state_file','fs2000_steadystate_initval');

%% compare the solutions
mwssf_a.print_solution
mwssf_b.print_solution

%% now we repeat the same exercise but with the steady state model inside the model file
% since we would like to do the two cases above in an elegant way, we
% introduce a switch called 'approximate'. You can call it what you want
% but it has to be issued while 'rising' the model. Not after.

% to continue in our elegance we initialize a vector of models
mv=rise.empty(0);
approximate=[0,1];
for approx=1:numel(approximate)
    approx_flag=struct('approximate',approximate(approx));
    % when approximate==0, we read the exact steady state solution
    % when approximate==1, we read an approximation of the steady state
    % solution. in both cases, the computed values will be used as an
    % initial guess for the computation of the steady state, unless we add
    % the attribute 'imposed' to 'steady_state_model', in which case rise
    % does not check that the initial guess actually solves the steady
    % state.
    mv(approx,1)=rise('fs2000_b','rise_flags',approx_flag);
end

%% we solve both models simultaneously
mv=solve(mv);

%% we print the solution
print_solution(mv)

%% a word of warning
% you might be tempted to or make the mistake of having both a
% steady_state_model block and a steady state file. RISE will use the
% steady state file and ignore the steady_state_model block!!! If you don't
% like this behavior, shoot an email to junior.maih AT gmail.com and let me
% know about your arguments for doing it differently.

author: jmaih


text:
function msm_to_mm_test10 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST10 tests MSM_TO_MM_ARRAY_REAL_SYMMETRIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST10\n' );
  fprintf ( 1, '  Convert an MSM to MM array real symmetric format.\n' );

  output_filename = 'msm_to_mm_test10.mm';

  a = r8mat_indicator ( 4, 4 );
  a = a + a';
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'array', 'real', 'symmetric' );

  return
end


text:
function [F1,F2]=framepair(ftype,g1,g2,varargin)
%FRAMEPAIR  Construct a new frame
%   Usage: [F1,F2]=framepair(ftype,g1,g2,...);
%
%   `[F1,F2]=framepair(ftype,g1,g2,...)` constructs two new frame objects
%   *F1* and *F2* of the same type *ftype* using the windows *g1* and *g2*.
%   The windows are specific to choosen frame type. See the help on *frame*
%   for the windows and arguments.
%
%   This function makes it easy to create a pair of canonical dual frames:
%   simply specify `'dual'` as window if one frame should be the dual of the
%   other.
%
%   This is most easily explained through some examples. The following
%   example creates a Gabor frame for real-valued signals with a Gaussian
%   analysis window and its canonical dual frame as the synthesis frame:::
%
%      f=greasy;
%      [Fa,Fs]=framepair('dgtreal','gauss','dual',20,294);
%      c=frana(Fa,f);
%      r=frsyn(Fs,c);
%      norm(f-r)
%
%   The following example creates a Wilson basis with a Gaussian
%   synthesis window, and its canonical dual frame as the analysis
%   frame::
%
%     [Fa,Fs]=framepair('dwilt','dual','gauss',20);
%
%   See also: frame, framedual, frametight

complainif_notenoughargs(nargin,3,'FRAMEPAIR');

ftype=lower(ftype);

if ~strcmp(g1,'dual')
    F1=frame(ftype,g1,varargin{:});
end;

if ~strcmp(g2,'dual')
    F2=frame(ftype,g2,varargin{:});
end;

if strcmp(g1,'dual')
    F1=framedual(F2);
end;

if strcmp(g2,'dual')
    F2=framedual(F1);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function t_most_30b_3_1_0(quiet)
%T_MOST_30B_3_1_0  Tests for MOST.

%   MOST
%   Copyright (c) 2009-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

if nargin < 1
    quiet = 0;
end

n_tests = 39;

t_begin(n_tests, quiet);

casename = 't_case30_most';
fudging = struct( ...       %% paramters for fudging reserve contract for sopf2
    'fudge',    0.05, ...   %% initial value (MW)
    'step',     0.01, ...   %% if necessary, increase by this amount and retry (MW)
    'lim',      0.1);       %% upper limit (MW), give up if no convergence
                            %% with fudge equal to this limit

%% options
mpopt = mpoption('verbose', 0, 'out.all', 0);
mpopt = mpoption(mpopt, 'opf.violation', 5e-7, 'mips.comptol', 5e-8);
if have_feature('linprog')
    if have_feature('linprog_ds')
        mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
    else
        mpopt = mpoption(mpopt, 'linprog.Algorithm', 'simplex');
    end
end
mpoptac = mpoption(mpopt, 'model', 'AC');
mpoptdc = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'most.solver', 'DEFAULT');

%% turn off warnings
s7 = warning('query', 'MATLAB:nearlySingularMatrix');
s6 = warning('query', 'MATLAB:nearlySingularMatrixUMFPACK');
warning('off', 'MATLAB:nearlySingularMatrix');
warning('off', 'MATLAB:nearlySingularMatrixUMFPACK');

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
    CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
    CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
    CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
    CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
    CT_MODCOST_X] = idx_ct;

%% reserve and delta offers
xgd_table.colnames = {
    'PositiveActiveReservePrice', ...
            'PositiveActiveReserveQuantity', ...
                    'NegativeActiveReservePrice', ...
                            'NegativeActiveReserveQuantity', ...
                                    'PositiveActiveDeltaPrice', ...
                                            'NegativeActiveDeltaPrice', ...
};
xgd_table.data = [
    10.1    15      10.0    15      0.1     0.0;
    10.3    30      10.2    30      0.3     0.2;
    10.5    20      10.4    20      0.5     0.4;
    10.7    25      10.6    25      0.7     0.6;
    20.1    25      20.0    25      60.1    60.0;
    20.3    15      20.2    15      15.1    15.0;
    20.5    30      20.4    30      60.3    60.2;
    20.7    15      20.6    15      15.3    15.2;
    30.1    15      30.0    15      60.3    60.4;
    30.3    30      30.2    30      30.1    30.0;
    30.5    25      30.4    25      60.7    60.6;
    30.7    30      30.6    30      30.3    30.2;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
];

%% contingency table
% label probty  type        row column      chgtype newvalue
contab = [
%     1   0.002   CT_TBRCH    1   BR_STATUS   CT_REP  0;      %% line 1-2
%     2   0.002   CT_TBRCH    2   BR_STATUS   CT_REP  0;      %% line 1-3, all power from gen1 flows via gen2
%     3   0.002   CT_TBRCH    3   BR_STATUS   CT_REP  0;      %% line 2-4, a path to loads @ buses 7 & 8
%     4   0.002   CT_TBRCH    5   BR_STATUS   CT_REP  0;      %% line 2-5, a path to loads @ buses 7 & 8
%     5   0.002   CT_TBRCH    6   BR_STATUS   CT_REP  0;      %% line 2-6, a path to loads @ buses 7 & 8
%     6   0.002   CT_TBRCH    36  BR_STATUS   CT_REP  0;      %% line 28-27, tie line between areas 1 & 3
%     7   0.002   CT_TBRCH    15  BR_STATUS   CT_REP  0;      %% line 4-12, tie line between areas 1 & 2
%     8   0.002   CT_TBRCH    12  BR_STATUS   CT_REP  0;      %% line 6-10, tie line between areas 1 & 3
%     9   0.002   CT_TBRCH    14  BR_STATUS   CT_REP  0;      %% line 9-10, tie line between areas 1 & 3
%     10  0.002   CT_TGEN     1   GEN_STATUS  CT_REP  0;      %% gen 1 at bus 1
%     11  0.002   CT_TGEN     2   GEN_STATUS  CT_REP  0;      %% gen 2 at bus 2
%     12  0.002   CT_TGEN     3   GEN_STATUS  CT_REP  0;      %% gen 3 at bus 22
%     13  0.002   CT_TGEN     4   GEN_STATUS  CT_REP  0;      %% gen 4 at bus 27
%     14  0.002   CT_TGEN     5   GEN_STATUS  CT_REP  0;      %% gen 5 at bus 23
%     15  0.002   CT_TGEN     6   GEN_STATUS  CT_REP  0;      %% gen 6 at bus 13
%     20  0.010   CT_TGEN     0   PMIN        CT_REL  1.1;    %% 10% load increase
%     20  0.010   CT_TGEN     0   QMIN        CT_REL  1.1;
%     21  0.010   CT_TGEN     0   PMIN        CT_REL  0.9;    %% 10% load decrease
%     21  0.010   CT_TGEN     0   QMIN        CT_REL  0.9;
];
clist = [];
nc = length(clist);

%% load the case
mpc = loadcase(casename);
gbus = mpc.gen(:, GEN_BUS);
mpc.gen(:, RAMP_30) = Inf;      %% avoid binding ramp reserves in period 1

%%-----  get OPF results  -----
rdc = rundcopf(mpc, mpoptdc);
% rac = runopf(mpc, mpoptac);
% save t_most4_soln rdc rac -v6
% s = load('t_most4_soln');
s.rdc = rdc;
% s.rac = rac;

%%-----  set up data for DC run (most)  -----
ng = size(mpc.gen, 1);      %% number of gens
nt = 3;
xgd = loadxgendata(xgd_table, mpc);
md = loadmd(mpc, nt, xgd);

%%-----  do DC run (most)  -----
r = most(md, mpopt);

%%-----  test the results  -----
t = 'success1';
t_ok(s.rdc.success, t);
t = 'success2';
t_is(r.QP.exitflag, 1, 12, t);

t = 'f';
t_is(r.results.f/sum(r.StepProb), s.rdc.f, 4, t);

for tt = 1:nt

t = sprintf('(t=%d) Pg : base', tt);
t_is(r.flow(tt,1,1).mpc.gen(:, PG), s.rdc.gen(:, PG), 5, t);

t = sprintf('(t=%d) gen : base', tt);
t_is(r.flow(tt,1,1).mpc.gen(:,1:MU_PMIN), s.rdc.gen(:,1:MU_PMIN), 3, t);

t = sprintf('(t=%d) energy prices', tt);
t_is(r.results.GenPrices(:,tt), s.rdc.bus(gbus, LAM_P), 6, t);

t = sprintf('(t=%d) Pc', tt);
t_is(r.results.Pc(:,tt), s.rdc.gen(:, PG), 4, t);

t = sprintf('(t=%d) Gmin', tt);
t_is(r.results.Pc(:,tt) - r.results.Rpm(:,tt), s.rdc.gen(:, PG), 4, t);

t = sprintf('(t=%d) Gmax', tt);
t_is(r.results.Pc(:,tt) + r.results.Rpp(:,tt), s.rdc.gen(:, PG), 4, t);

t = sprintf('(t=%d) upward contingency reserve quantities', tt);
t_is(r.results.Rpp(:,tt), zeros(ng, 1), 4, t);

t = sprintf('(t=%d) downward contingency reserve quantities', tt);
t_is(r.results.Rpm(:,tt), zeros(ng, 1), 4, t);

t = sprintf('(t=%d) upward contingency reserve prices', tt);
t_is(r.results.RpmPrices(:,tt), xgd.NegativeActiveReservePrice, 6, t);

t = sprintf('(t=%d) downward contingency reserve prices', tt);
t_is(r.results.RpmPrices(:,tt), xgd.NegativeActiveReservePrice, 6, t);

t = sprintf('(t=%d) Rpmax_pos', tt);
vv = r.om.get_idx();
Rpmax_pos = (r.QP.lambda.upper(vv.i1.Rpp(tt):vv.iN.Rpp(tt)) - r.QP.lambda.lower(vv.i1.Rpp(tt):vv.iN.Rpp(tt))) / mpc.baseMVA;
t_is(Rpmax_pos, zeros(ng, 1), 6, t);

t = sprintf('(t=%d) Rpmax_neg', tt);
Rpmax_neg = (r.QP.lambda.upper(vv.i1.Rpm(tt):vv.iN.Rpm(tt)) - r.QP.lambda.lower(vv.i1.Rpm(tt):vv.iN.Rpm(tt))) / mpc.baseMVA;
t_is(Rpmax_neg, zeros(ng, 1), 6, t);

end

% g1 = s.rdc.base.gen(:, PG);
% g2 = r.flow(1,1,1).mpc.gen(:, PG);
% for k = 1:nc
%     g1 = [ g1 s.rdc.cont(k).gen(:, PG) ];
%     g2 = [ g2 r.flow(1,1,k+1).mpc.gen(:, PG) ];
% end
% [m,n] = size(g1);
% for j = 1:n
%     fprintf('\n');
%     for i = 1:m
%         fprintf('%9.2f  %9.2f\n', g1(i,j), g2(i,j));
%     end
% end

%%-----  do AC run (most)  -----
%mostac;

%% turn warnings back on
warning(s7.state, 'MATLAB:nearlySingularMatrix');
warning(s6.state, 'MATLAB:nearlySingularMatrixUMFPACK');

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function stroud_test45 ( )

%*****************************************************************************80
%
%% TEST45 tests TORUS_1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_3D_INDEX;

  r1 = 0.5;
  r2 = 1.0;
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST45\n' );
  fprintf ( 1, '  TORUS_1 approximates integrals on a torus.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The degree N will be varied.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Inner radius = %f\n', r1 );
  fprintf ( 1, '  Outer radius = %f\n', r2 );
  fprintf ( 1, '  Area = %f\n', torus_area_3d ( r1, r2 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    F(X)  ' );
  for j = 0 : 2 : 8
    fprintf ( 1, '  %6d      ', 2^j );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );

  num = function_3d_num ( );

  for i = 1 : num

    FUNC_3D_INDEX = i;

    for j = 1 : 5

      j2 = 2 * ( j - 1 );
      n = 2^j2;
      result(j) = torus_1 ( 'function_3d', r1, r2, n );

    end

    fname = function_3d_name ( i );
    fprintf ( 1, '  %s', fname );
    for j = 1 : 5
      fprintf ( 1, '  %12f', result(j) );
    end
    fprintf ( 1, '\n' );

  end

  return
end


text:
%GETUSER Return the user field of a mapping
%
%    USERFIELD = GETUSER(W,FIELD)
%
% INPUT
%   W          Mapping
%   FIELD      Character string: name of structure field of USER field of W.
%
% OUTPUT
%   USERFIELD  Requested field
%
% If FIELD is is the empty string ('') the entire user field is returned.
%
% If the requested field does not exist USERFIELD = [];
%
% Note the the USER field of mappings was originally intended for a user
% defined description of mappings. Later its usage was extended to a field
% for storing general information on mappings. For that reason 'old'
% mappings without a structure in the user field are transformed such that
% this information is stored in a subfield USER in the user field. It can
% be retrieved by GETUSER(W).

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%% Copyright (C) 2014-2017, 2019 Colin B. Macdonald
%%
%% 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
%% @defop  Method   @@sym {horzcat} {(@var{x}, @var{y}, @dots{})}
%% @defopx Operator @@sym {[@var{x}, @var{y}, @dots{}]} {}
%% @defopx Operator @@sym {[@var{x} @var{y} @dots{}]} {}
%% Horizontally concatentate symbolic arrays.
%%
%% Example:
%% @example
%% @group
%% A = sym([1 2; 3 4])
%%   @result{} A = (sym 2×2 matrix)
%%       ⎡1  2⎤
%%       ⎢    ⎥
%%       ⎣3  4⎦
%%
%% [A A A]
%%   @result{} (sym 2×6 matrix)
%%       ⎡1  2  1  2  1  2⎤
%%       ⎢                ⎥
%%       ⎣3  4  3  4  3  4⎦
%% @end group
%% @end example
%% @seealso{@@sym/vertcat, @@sym/cat}
%% @end defop

function h = horzcat(varargin)

  % special case for 0x0 but other empties should be checked for
  % compatibilty
  cmd = {
          '_proc = []'
          'for i in _ins:'
          '    if i is None or not i.is_Matrix:'
          '        _proc.append(sp.Matrix([[i]]))'
          '    else:'
          '        if i.shape == (0, 0):'
          '            pass'
          '        else:'
          '            _proc.append(i)'
          'return sp.MatrixBase.hstack(*_proc),'
          };

  for i = 1:nargin
    varargin{i} = sym(varargin{i});
  end
  h = pycall_sympy__ (cmd, varargin{:});

end

%!test
%! % basic
%! syms x
%! A = [x x];
%! B = horzcat(x, x);
%! C = horzcat(x, x, x);
%! assert (isa (A, 'sym'))
%! assert (isa (B, 'sym'))
%! assert (isa (C, 'sym'))
%! assert (isequal (size(A), [1 2]))
%! assert (isequal (size(B), [1 2]))
%! assert (isequal (size(C), [1 3]))

%!test
%! % basic, part 2
%! syms x
%! A = [x 1];
%! B = [1 x];
%! C = [1 2 x];
%! assert (isa (A, 'sym'))
%! assert (isa (B, 'sym'))
%! assert (isa (C, 'sym'))
%! assert (isequal (size(A), [1 2]))
%! assert (isequal (size(B), [1 2]))
%! assert (isequal (size(C), [1 3]))

%!test
%! % row vectors
%! a = [sym(1) 2];
%! b = [sym(3) 4];
%! assert (isequal ( [a b] , [1 2 3 4]  ))
%! assert (isequal ( [a 3 4] , [1 2 3 4]  ))
%! assert (isequal ( [3 4 a] , [3 4 1 2]  ))
%! assert (isequal ( [a [3 4]] , [1 2 3 4]  ))
%! assert (isequal ( [a sym(3) 4] , [1 2 3 4]  ))
%! assert (isequal ( [a [sym(3) 4]] , [1 2 3 4]  ))

%!test
%! % col vectors
%! a = [sym(1); 2];
%! b = [sym(3); 4];
%! assert (isequal ( [a b] , [1 3; 2 4]  ))
%! assert (isequal ( [a b a] , [1 3 1; 2 4 2]  ))

%!test
%! % empty vectors
%! v = sym(1);
%! a = [v []];
%! assert (isequal (a, v))
%! a = [[] v []];
%! assert (isequal (a, v))
%! a = [v [] []];
%! assert (isequal (a, v))

%!test
%! % more empty vectors
%! v = [sym(1) sym(2)];
%! q = sym(ones(1, 0));
%! assert (isequal ([v q], v))

%!error <ShapeError>
%! v = [sym(1) sym(2)];
%! q = sym(ones(3, 0));
%! w = horzcat(v, q);

%!test
%! % issue #700
%! A = sym ([1 2]);
%! B = simplify (A);
%! assert (isequal ([B A], [A B]))

repo: octsympy


text:
function DEM = spm_MDP_DEM(DEM,demi,O,o)
% auxiliary (link) function for mixed hierarchical (MDP/DEM) models
% FORMAT DEM = spm_MDP_DEM(DEM,demi,O,o)
%
% DEM      - DEM structure
% demi     - mapping from discrete outcomes to hidden causes
%  demi.C  - cell array of true causes for each combination of outcomes
%            the appropriate array is then placed in DEM.C
%  demi.U  - cell array of hidden causes for each combination of outcomes
%            the Bayesian model average is placed in DEM.U
% O{g}     - cell array of priors over discrete outcomes
% o(g x 1) - vector of true outcomes
%
% completes the following fields:
%   DEM.X{g} - posterior probability over g models and t times
%
% This routine performs a Bayesian model comparison using (DEM) Bayesian
% filtering and places the results in fields of the DEM structure; so that
% MDP schemes can pick them up as likelihood terms in the next hierarchical
% level. The outcomes of the (discrete) MDP scheme at the superordinate
% level specify the hidden causes at the current level. These enter as
% Bayesian model averages of the continuous causes. The resulting
% posterior over hidden causes then furnishes the posterior over outcomes
% using Bayesian model reduction, based on the free energy accumulated
% (integrated) over time. This free energy is supplemented with the prior
% over discrete outcomes; thereby constituting a posterior over outcomes.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_DEM.m 7560 2019-03-29 14:39:38Z thomas $

% evaluate true values and priors over causes given discrete states
%--------------------------------------------------------------------------
po    = spm_cross(O);
ind   = num2cell(o);
DEM.C = demi.C{ind{:}};
DEM.U = 0;
for i = 1:size(po,1)
    for j = 1:size(po,2)
        for k = 1:size(po,3)
            DEM.U = DEM.U + demi.U{i,j,k}*po(i,j,k);
        end
    end
end

% integrate system to generate data
%--------------------------------------------------------------------------
DEM.db = 0;
DEM    = spm_ADEM(DEM);

% posterior probability over discrete states models
%==========================================================================
nv    = size(DEM.U,1);
nt    = size(DEM.U,2);
ic    = (1:nv) + size(DEM.qU.C{1},1) - nv;
P     = spm_zeros(po);
for i = 1:ndims(P)
    x{i} = ones(size(P,i),1);
end

% prior potential
%--------------------------------------------------------------------------
F     = log(po);

% and accumulate log evidence
%--------------------------------------------------------------------------
for t = 1:nt

    % posteriors and priors for this time point
    %----------------------------------------------------------------------
    qE    = DEM.qU.v{end}(:,t);
    qC    = DEM.qU.C{t}(ic,ic);
    pE    = DEM.U(:,t);
    pC    = DEM.pU.C{t}(ic,ic);
    gt    = t > 4;
    for i = 1:size(po,1)
        for j = 1:size(po,2)
            for k = 1:size(po,3)
                if po(i,j,k) < exp(-3)
                    F(i,j,k) = - exp(64);
                else
                    rE       = demi.U{i,j,k}(:,t);
                    F(i,j,k) = F(i,j,k) + gt*spm_log_evidence(qE,qC,pE,pC,rE,pC);
                end
            end
        end
    end

    % Bayesian model comparison
    %----------------------------------------------------------------------
    P(:)  = spm_softmax(F(:));

    % marginal posteriors
    %----------------------------------------------------------------------
    for g = 1:numel(O)
        X{g}(:,t) = spm_dot(P,x,g);
    end
end

% return  probability over models (i.e., outcomes at subordinate level)
%--------------------------------------------------------------------------
DEM.X = X;

author: spm


text:
function varargout = mesh(varargin)
%MESH   3-D mesh surface of a SPHEREFUN.
%   MESH is not supported for SPHEREFUN objects, and so returns an error.
%
% See also PLOT, SURF.

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

[varargout{1:nargout}] = mesh@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [A, b, Aeq, beq, lb, ub] = snd_l_constraints(k, net_load_data, varargin)

    u0_ref = cell2mat(varargin(2));
    f = 1.2;
    g = 0.8;
    p = 1.4;
    A   = []; %A*u<=b
    b   = [];
    Aeq = [1 1 1]; %A*u = b
    beq = net_load_data(k);

    if k ~= 12
        if u0_ref(1,k) * u0_ref(1,k+1) <= 0
            if u0_ref(1,k)>0
                ub = [u0_ref(1,k)*f];
                lb = [u0_ref(1,k+1)*f];
            else
                ub = [u0_ref(1,k+1)*f];
                lb = [u0_ref(1,k)*f];
            end
        elseif u0_ref(1,k)>0 && u0_ref(1,k+1)>0
            if u0_ref(1,k)>u0_ref(1,k+1)
                ub = [u0_ref(1,k)*f];
                lb = [u0_ref(1,k+1)*g];
            else
                ub = [u0_ref(1,k+1)*f];
                lb = [u0_ref(1,k)*g];
            end
        elseif u0_ref(1,k)<0 && u0_ref(1,k+1)<0
            if u0_ref(1,k)>u0_ref(1,k+1)
                ub = [u0_ref(1,k)*g];
                lb = [u0_ref(1,k+1)*f];
            else
                ub = [u0_ref(1,k+1)*g];
                lb = [u0_ref(1,k)*f];
            end
        end
    else
        if u0_ref(1,k) >=0
            ub = [u0_ref(1,k)*f];
            lb = [u0_ref(1,k)*g];
        else
            ub = [u0_ref(1,k)*g];
            lb = [u0_ref(1,k)*f];
        end

    end

    if u0_ref(2,k)>0 %
        ub  = [ub, u0_ref(2,k)*p,  10];
        lb  = [lb, 0, -10];
    else
        ub  = [ub, 0,  10];
        lb  = [lb, u0_ref(2,k)*p, -10];
    end

end

author: juchengquan


text:
function DCM = spm_meta_model(DCM)
% Meta-modelling of Bayes-optimal responses (Newton's method)
% FORMAT DCM = spm_meta_model(DCM)
%
% store estimates in DCM
%--------------------------------------------------------------------------
% DCM.M   - meta-model specification
%      M: [1 x m struct]   - hierarchical inference model (cf DEM.M)
%      G: [1 x s struct]   - generative process (for spm_ADEM or spm_ALAP)
%      U: [n x N double]   - n prior beliefs over N samples
%     pE: [1 x 1 struct]   - prior expectation of meta-model parameters
%     pC: [1 x 1 struct]   - prior variance of meta-model parameters
%
% DCM.xY  - data structure
%      y: [N x p double]   - N samples of a p-variate response
%     X0: [N x q double]   - q-variate confounds
%     dt: [1 x 1 double]   - size of time bin for each sample
%      Q: {[N x N double]} - precision component[s]
%
% DCM.xU  - input structure
%      u: [r x N double]   - r-variate input (hidden causes G in DEM)
%
% Computes (and stores in DCM_MM_???)
%--------------------------------------------------------------------------
% DCM.DEM - Inference (with MAP parameters)
% DCM.Ep  - conditional expectation
% DCM.Cp  - conditional covariances
% DCM.Eh  - conditional log-precision
% DCM.Ey  - conditional response
% DCM.F   - log-evidence
%
% This routine illustrates Bayesian meta modelling - the Bayesian inversion
% of a model of a Bayesian observer. This requires the specification of two
% models: An inference model used by the subject (specified by a DEM
% structure) and a meta-model (specified by a DCM structure). The inference
% model is completed by a response model to furnish the meta-model; where
% the response model takes the output of the (active) inference scheme
% specified by the DEM and generates an observed (behavioural or
% neurophysiological) response. Crucially either the inference model or
% the response model or both can have free parameters - that are optimised
% using Bayesian nonlinear system identification in the usual way.
%
% Although this routine is a function, it is expected that people will fill
% in the model-specific parts in a local copy, before running it.  The
% current example uses a model of slow pursuit and generates synthetic data
% (responses) to illustrate how it works. To replace these simulated data
% with real data, simply specify the DCM.xY (and xU fields) with
% empirical values. If other fields do not exist, exemplar fields will be filled in.
%
% The conditional density of the parameters and F values (log-evidence) can
% be used in the usual way for inference on parameters or Bayesian model
% comparison (as for other DCMs)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_meta_model.m 5788 2013-12-06 20:08:57Z karl $

% INFERENCE MODEL: DEM.M and DEM.G
%==========================================================================
try
    M = DCM.M;

catch

    % create illustrative (M,G) if not specified (EXAMPLE)
    %----------------------------------------------------------------------

    % hidden causes and states
    %======================================================================
    % x    - hidden states:
    %   x.o(1) - oculomotor angle
    %   x.o(2) - oculomotor velocity
    %   x.x(1) - target angle - extrinsic coordinates
    %
    % v    - causal states: force on target
    %
    % g    - sensations:
    %   g(1) - oculomotor angle (proprioception)
    %   g(2) - oculomotor velocity
    %   g(:) - visual input - intrinsic coordinates
    %----------------------------------------------------------------------

    % Set-up
    %======================================================================
    M(1).E.s = 1/2;                               % smoothness
    M(1).E.n = 4;                                 % order of
    M(1).E.d = 1;                                 % generalised motion

    % angular frequency of target motion (per time bin)
    %----------------------------------------------------------------------
    w  = 2*pi/32;

    % occlusion and visual mapping functions
    %----------------------------------------------------------------------
    occ = @(x) x < 1/2;
    vis = @(x) exp(-((-8:8)' - x.x(1) + x.o(1)).^2)*occ(x.x);

    % Generative model (M) (sinusoidal movement)
    %======================================================================
    % Endow the model with internal dynamics (a simple oscillator) so that is
    % recognises and remembers the trajectory to anticipate jumps in rectified
    % sinusoidal motion.

    % slow pursuit following with (second order) generative model
    %----------------------------------------------------------------------
    x.o = [0;0];                                  % motor angle & velocity
    x.x = 0;                                      % target location

    % level 1: Displacement dynamics and mapping to sensory/proprioception
    %----------------------------------------------------------------------
    M(1).f = @(x,v,P) [x.o(2); (v - x.o(1))/4 - x.o(2)/2; v - x.x];
    M(1).g = @(x,v,P) [x.o; vis(x)];
    M(1).x = x;                                   % hidden states
    M(1).V = exp(4);                              % error precision
    M(1).W = exp(4);                              % error precision

    % level 2: With hidden (memory) states
    %----------------------------------------------------------------------
    M(2).f  = @(x,v,P)[x(2); -x(1)]*v;
    M(2).g  = @(x,v,P) x(1);
    M(2).x  = [0; 0];                             % hidden states
    M(2).V  = exp(4);                             % error precision
    M(2).W  = exp(4);                             % error precision

    % level 3: Encoding frequency of memory states (U)
    %----------------------------------------------------------------------
    M(3).v = 0;
    M(3).V = exp(4);

    % generative model (G) (with no noise)
    %======================================================================

    % first level
    %----------------------------------------------------------------------
    G(1).f = @(x,v,a,P) [x.o(2); a/4 - x.o(2)/8; v - x.x];
    G(1).g = @(x,v,a,P) [x.o; vis(x)];
    G(1).x = x;                                  % hidden states
    G(1).U = sparse(1,[1 2],[1 1],1,19)*exp(4);  % motor gain

    % second level
    %----------------------------------------------------------------------
    G(2).v = 0;                                  % exogenous force
    G(2).a = 0;                                  % action force

end

% Check generative model
%--------------------------------------------------------------------------
M      = spm_DEM_M_set(M);

% Experimental input (C)
%==========================================================================
try
    xU    = DCM.xU;
    N     = size(xU.u,2);

catch

    % Sine wave cause (EXAMPLE)
    %----------------------------------------------------------------------
    N     = 64;                                 % length of data sequence
    xU.dt = 16;                                 % time step (ms)
    xU.u  = sin((1:N)*w).*((1:N) > 16);         % sinusoidal target motion
end

% Priors (U)
%==========================================================================
try
    U = DCM.U;

catch

    % (EXAMPLE) prior beliefs about target frequency (w)
    %----------------------------------------------------------------------
    U = zeros(1,N) + w;
end

% Data and confounds
%==========================================================================
try
    xY = DCM.xY;

catch

    % Generate simulated data (EXAMPLE)
    %----------------------------------------------------------------------

    % meta-model and TRUE parameters
    %----------------------------------------------------------------------
    MM.M    = M;
    MM.G    = G;
    MM.U    = U;

    P.G     = 2;
    P.W     = 4;
    [y DEM] = IS(P,MM,xU);

    % and simulate observation noise
    %----------------------------------------------------------------------
    xY.y    = y + spm_Q(1/2,N,1)*randn(N,1)/16;

    spm_figure('GetWin','Simulated respsone');
    spm_DEM_qU(DEM.qU,DEM.pU)

end

% DCT confounds
%==========================================================================
xY.X0  = spm_dctmtx(N,1);

% and [serial] correlations (precision components) AR model
%--------------------------------------------------------------------------
xY.Q   = {spm_Q(1/2,N,1)};

% META-MODEL parameters
%==========================================================================
% These parameterise the inference scheme function (IS) at the bottom of
% this script - this functions defines how the parameters are used and
% therefore optimised. Change this function to specify the meta-model

try
    pE    = DCM.pE;
    pC    = DCM.pC;

catch

    % prior expectations (EXAMPLE)
    %----------------------------------------------------------------------
    pE.G  = 1;
    pE.W  = 2;

    % prior covariance (EXAMPLE)
    %----------------------------------------------------------------------
    pC.G  = 1;
    pC.W  = 1;

end

% This completes the meta-model specification. The fields are now assembled
% and passed to spm_nlsi_GN (nonlinear system identification using Gauss-
% Newton-like gradient ascent).
%==========================================================================

% META-MODEL (MM)
%==========================================================================
MM.M   = M;
MM.G   = G;
MM.U   = U;

% hyperpriors (assuming about 99% signal to noise)
%--------------------------------------------------------------------------
hE     = 6 - log(var(spm_vec(xY.y)));
hC     = exp(-4);

% Meta-model
%--------------------------------------------------------------------------
MM.IS  = @IS;
MM.pE  = pE;
MM.pC  = pC;
MM.hE  = hE;
MM.hC  = hC;

% model inversion
%==========================================================================
MM.Nmax      = 16;
[Ep,Cp,Eh,F] = spm_nlsi_GN(MM,xU,xY);

% integrate (A)DEM scheme with MAP etimates
%--------------------------------------------------------------------------
[Ey,DEM] = IS(Ep,MM,xU);

% store estimates in DCM
%--------------------------------------------------------------------------
DCM.DEM = DEM;                  % Inference scheme (with MAP parameters)
DCM.M   = MM;                   % meta-model
DCM.xY  = xY;                   % data structure
DCM.xU  = xU;                   % input structure
DCM.Ep  = Ep;                   % conditional expectation
DCM.Cp  = Cp;                   % conditional covariances
DCM.Eh  = Eh;                   % conditional log-precision
DCM.Ey  = Ey;                   % conditional response
DCM.F   = F;                    % log-evidence

% save
%--------------------------------------------------------------------------
save(sprintf('DCM_MM_%s',date),'DCM', spm_get_defaults('mat.format'));

% and plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Parameter estimates');
subplot(2,1,1)
spm_plot_ci(Ep,Cp), hold on
Pp       = spm_vec(P);
pE       = spm_vec(pE);
plot(Pp(1),Pp(2),'.r' ,'MarkerSize',32), hold on
plot(pE(1),pE(2),'.g','MarkerSize',32), hold off
axis([0 4 0 8])
xlabel('Gain','FontSize',12)
ylabel('Precision','FontSize',12)
title('Posterior (blue), prior (green) and true (red) value','FontSize',16)

return

% inference scheme x = IS(p,M,U)
%==========================================================================
function [y,DEM] = IS(P,M,U)

% parameterise inference model (EXAMPLE)
%--------------------------------------------------------------------------
M.M(2).W = exp(P.W);

% reset random number generator (for action schemes)
%--------------------------------------------------------------------------
rng('default')

% Bayesian inference (spm_DEM, spm_LAP, spm_ADEM or spm_ALAP)
%--------------------------------------------------------------------------
DEM.M  = M.M;
DEM.G  = M.G;
DEM.U  = M.U;
DEM.C  = U.u;
DEM    = spm_ADEM(DEM);

% (Hidden) behavioural or electrophysiological response (EXAMPLE)
%--------------------------------------------------------------------------
x      = DEM.pU.x{1}(1,:);
y      = (P.G*x)';

return


text:
clear
close all

plotx2east

load martensite_single_grain

[grains,ebsd('indexed').grainId] =...
    calcGrains(ebsd('indexed'),'angle',3*degree);

cond_grains = ismember(grains.id,unique(ebsd(ids_of_interest).grainId));

%%
job = parentGrainReconstructor(ebsd,grains) % define job
job.useBoundaryOrientations = true;
clearvars grains ebsd

job.p2c = orientation.NishiyamaWassermann(job.p2c.CS,job.p2c.SS);

%%
job.calcHyperGraph2('threshold',5*degree,'c2c')

numIters = [3 10 50];
grainId = 1844;

for k = 1:length(numIters);
    if k == 1
        job.clusterHyperGraph2('numIter',numIters(k),'inflationPower',1,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    else
        job.clusterHyperGraph2('numIter',numIters(k) - numIters(k-1),'inflationPower',1,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    end
end

for k = 1:length(numIters);
    if k == 1
        job.clusterHyperGraph2('numIter',numIters(k),'inflationPower',1.05,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    else
        job.clusterHyperGraph2('numIter',numIters(k) - numIters(k-1),'inflationPower',1.05,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    end
end

%%
function graphplotfunc(job,cond_grains,grainId,k,numIters)

pOri = variants(job.p2c, job.grains.meanOrientation, job.votes.parentId(:,1));

nextAxis
plot(job.grains(cond_grains),...
    pOri(cond_grains))
lims = job.ebsd.extent;
xlim([lims(1) lims(2)])
ylim([lims(3) lims(4)])

text(lims(1)+5,lims(4)-10, [num2str(numIters(k)) ' iterations'],...
    'HorizontalAlignment','left',...
    'FontSize',20, 'FontWeight','bold','BackgroundColor', 'white')

end

author: mtex-toolbox

repo: mtex

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


text: % VL_IHASHSUM

author: jianxiongxiao

repo: ProfXkit


text:
function Y = tapas_uniqc_create_shepp_logan_4d()
% Creates 4D Shepp-Logan Phantom for tapas_uniqc_slider4d demo
%
%   output = tapas_uniqc_create_shepp_logan_4d(input)
%
% IN
%
% OUT
%
% EXAMPLE
%   tapas_uniqc_create_shepp_logan_4d
%
%   See also

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

% parameter specs, create shepp logan phantom
nX = 64;
nY = nX;
nSli = 32;
nDyn = 20;

P = phantom('Modified Shepp-Logan',nX);

% create different slices via shift
Y = zeros(nX,nY,nSli,1);
for iSli = 1:nSli
    Y(:,:,iSli,1) = circshift(P, iSli*floor(nX/nSli));
end

% replicate over number of dynamics and add some noise to make them
% different
Y = repmat(Y, [1 1 1 nDyn]) + 0.05*max(Y(:))*randn(nX,nY,nSli,nDyn);

author: translationalneuromodeling

repo: tapas


text:
function sigma = IRnormest(A, b)
%  We sometimes need an estimate of the 2-norm of the matrix A.  We want to
%  allow for user-defined objects and function handles that implement matrix
%  vector multiplication with A.  So we cannot use MATLAB's built-in
%  normest, nor can we use svds.  So we'll use a few iterations of our
%  Lanczos bidiagonalization as implemented in HyBR to get an estimate of
%  the largest singular value.
% HyBRoptions = HyBRset;
% HyBRoptions = HyBRset(HyBRoptions, 'InSolv', 'Tikhonov', 'RegPar', 0, 'Iter', 5, 'Reorth', 'on', 'verbosity', 'off');
% [~, HyBRout] = HyBR_for_IRtools(A, b, [], HyBRoptions, 'off', 'off');
% sigma = max(svd(HyBRout.B));

optnrm.RegParam = 0;
optnrm.MaxIter = 5;
optnrm.Reorth = 'on';
optnrm.DecompOut = 'on';
optnrm.IterBar = 'off';
optnrm.verbosity = 'off';
[~, infonrm] = IRhybrid_lsqr(A, b, optnrm);
sigma = max(svd(infonrm.B));

repo: IRtools


text:
function sMap = rsom_lininit(Distances, msize, varargin)

%RSOM_LININIT computes linear initialization for the RSOM.
%
% sM = rsom_lininit(D, msize, [argID, value, ...])
%
%  sM = rsom_lininit(D, msize, 'lattice', 'hexa');
%
% Input and output arguments:
%   D          (matrix) dissimilarity data for training, size nData x nData
%   msize      (vector) map size
%   [argID,    (string) See below.
%    value]    (varies)
%
%   sM         (struct) RSOM map struct, the initialized map
%
% Here are the valid argument IDs and corresponding values.
%   'lattice'     (string) map lattice: 'hexa' or 'rect'
%   'shape'       (string) map shape: 'sheet', 'cyl' or 'toroid'
%   'name'        (string) name of the RSOM struct
%
% For more help, try 'type rsom_lininit' or check out online documentation.
% See also RSOM_RANDINIT, RSOM_BATCHTRAIN.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% rsom_lininit
%
% PURPOSE
%
% Initializes the RSOM in the input space linearly.
%
% SYNTAX
%
%  sM = rsom_lininit(D, msize);
%  sM = rsom_lininit(...,'argID',value,...);
%
% DESCRIPTION
%
% The grid in the input space is initialized along the eigenvectors
% corresponding to the largest eigenvalues. These are calculated by
% classical MDS from the dissimilarity matrix D. If D is generated from
% euclidean data, squared distances should be provided i.e.
% (D)_ij = ||x_i - x_j||^2.
%
% REFERENCES
%
%   Barbara Hammer, Alexander Hasenfuss: Topographic Mapping of Large
%   Dissimilarity Data Sets. Neural Computation 22(9): 2229-2284 (2010)
%
% REQUIRED INPUT ARGUMENTS
%
%   D          (matrix) dissimilarity data for training, size nData x nData
%   msize      (vector) map size
%
% OPTIONAL INPUT ARGUMENTS
%
%  argID (string) Argument identifier string (see below).
%  value (varies) Value for the argument (see below).
%
%  The optional arguments can be given as 'argID',value -pairs.
%  The valid IDs and corresponding values are listed below.
%
%  Below is the list of valid arguments:
%   'lattice'     (string) map lattice: 'hexa' or 'rect'
%   'shape'       (string) map shape: 'sheet', 'cyl' or 'toroid'
%   'name'        (string) name of the RSOM struct
%
% EXAMPLES
%
%   sM = rsom_lininit(D, [10 10]);
%   sM = rsom_lininit(D, [10 10], 'lattice', 'hexa');
%
% SEE ALSO
%
%   rsom_randinit    Initialize a RSOM randomly
%   rsom_batchtrain  Train a RSOM

% Contributed to SOM Toolbox vs2, December 7th, 2012 by Alexander Schulz
% Copyright (c) Alexander Schulz
% http://www.cis.hut.fi/projects/somtoolbox/

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

%% Prase input
p = inputParser;
p.addRequired('Distances', @isnumeric);
p.addRequired('msize', @isnumeric);

p.addParamValue('lattice', 'rect', @ischar);
p.addParamValue('shape', 'sheet', @ischar);
p.addParamValue('name', '', @ischar);

p.parse(Distances, msize, varargin{:});

lattice = p.Results.lattice;
shape   = p.Results.shape;
name    = p.Results.name;
clear p;
mdim     = length(msize);
nData    = size(Distances,1);

%% Init and normalize
% create the topology struct
sTopol = som_topol_struct('msize', msize, ...
                          'lattice', lattice, ...
                          'shape', shape);

Coords = som_unit_coords(msize,'rect','sheet');
% normalize the neuron positions to mean 0
Coords = bsxfun(@minus, Coords, mean(Coords,1));

% compute mds embedding of the distance matrix
[V, e]=cmdscale(Distances.^(1/2));
% keep only the first mdim dimensions
V = V(:, 1:mdim);
e = e(1:mdim);
scale = std(V);

% scale Coords to the scale of the data
%Coords = bsxfun(@rdivide, Coords, std(Coords));
%Coords = bsxfun(@times, Coords, scale);
Coords = bsxfun(@rdivide, Coords, std(Coords)+1e-10);

% project neurons onto the data
% e is n * variance
standardDev = sqrt(e/nData);
V = bsxfun(@rdivide, V, nData*standardDev'+1e-10);
%V = bsxfun(@rdivide, V, sqrt(sum(V.^2, 2))+1e-10);

cCodebook = Coords * V';
cCodebook = cCodebook + 1/nData;

%% create the resulting struct
sTrain = som_train_struct('algorithm','lininit');

sMap = struct('type', 'rsom_map', ...
              'cCodebook', cCodebook, ...
              'topol', sTopol, ...
              'labels', cell(1), ...
              'neigh', 'gaussian', ...
              'trainhist', cell(1), ...
              'name', name);

sTrain = som_set(sTrain,'time',datestr(now,0));
sMap.trainhist = sTrain;

author: ilarinieminen

repo: SOM-Toolbox


text:
% SP_NURBS: Constructor of the class of a tensor-product spaces of NURBS.
%
%     sp = sp_nurbs (nurbs, msh)
%     sp = sp_nurbs (knots, degree, weights, msh)
%
% INPUTS:
%
%     nurbs:     nurbs structure from the NURBS toolbox (see nrbmak)
%     msh:       msh object that defines the quadrature rule (see msh_cartesian)
%     knots:     open knot vector (cell array of size [1, ndim])
%     degree:    nurbs polynomial degree (vector of size [1, ndim])
%     weights:   weights associated to the basis functions
%
% OUTPUT:
%
%    sp: object of the class sp_scalar, representing the discrete function space. See sp_scalar for details.
%
% Copyright (C) 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function sp = sp_nurbs (varargin)

  if (nargin <= 3)
    nurbs = varargin{1};
    msh   = varargin{2};
    if (nargin == 3)
      transform = varargin{3};
    else
      transform = 'grad-preserving';
    end

    knots   = nurbs.knots;
    degree  = nurbs.order - 1;
    weights = squeeze (nurbs.coefs(4, :, :, :));
  else
    knots   = varargin{1};
    degree  = varargin{2};
    weights = varargin{3};
    msh     = varargin{4};
    if (nargin == 5)
      transform = varargin{5};
    else
      transform = 'grad-preserving';
    end
  end

  sp = sp_scalar (knots, degree, weights, msh, transform);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%%***********************************************************************
%% nzlist: find the combined list of non-zero elements
%%         of Aj, j = 1:k, for each k,
%%         assuming that the Aj's are permuted such that
%%         A1 has the fewest nonzero elements, followed by A2, and so on.
%%
%% [isspA,nzlistA,nzlistAsum,isspAy,nzlistAy] = nzlist(blk,At,par)
%%
%%  isspA(p,k) = 1 if Apk is sparse, 0 if it is dense.
%%  nzlistA = px2 cell array.
%%            nzlistA{p,1}(k) is the starting row index (in C convention)
%%            in the 2-column matrix nzlistA{p,2} that
%%            stores the row and column index of the nonzero elements
%%            of Apk.
%%            nzlistA{p,1}(k) = inf if nnz(Apk) exceeds given threshold.
%%  nzlistAsum = px2 cell array.
%%            nzlistA{p,1}(k) is the starting row index (in C convention)
%%            in the 2-column matrix nzlistA{p,2} that
%%            stores the row and column index of the nonzero elements
%%            of Apk that are not already present
%%            in the combined list from Ap1+...Ap,k-1.
%%  nzlistAy = px1 cell array.
%%            nzlistAy{p} is a 2-column matrix that stores the
%%            row and column index of the nonzero elements of
%%            Ap,1+.... Ap,m.
%%            nzlistAy{p} = inf  if the number of nonzero elements
%%            exceeds a given threshold.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%***********************************************************************

function [isspA,nzlistA,nzlistAsum,isspAy,nzlistAy] = nzlist(blk,At,par)

spdensity   = par.spdensity;
smallblkdim = par.smallblkdim;
m = par.numcolAt;
%%
numblk = size(blk,1);
isspA = zeros(numblk,m);
nzlistA = cell(numblk,2);  nzlistAsum = cell(numblk,2);
isspAy = zeros(numblk,1);  nzlistAy = cell(numblk,1);
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if strcmp(pblk{1},'s') && ((max(pblk{2}) > smallblkdim) || (length(pblk{2}) <= 10))
        numblk = length(pblk{2});
        n = sum(pblk{2});
        n2 = sum(pblk{2}.*pblk{2});
        if (numblk == 1)
            nztol = spdensity*n;
            nztol2 = spdensity*n2/2;
        else
            nztol = spdensity*n/2;
            nztol2 = spdensity*n2/4;
        end
        nzlist1 = zeros(1,m+1); nzlist2 = [];
        nzlist3 = zeros(1,m+1); nzlist4 = []; breakyes = zeros(1,2);
        Asum = sparse(n,n);
        %%
        m1 = size(At{p,1},2);
        for k = 1:m1
            Ak = mexsmat(blk,At,1,p,k);
            nnzAk = nnz(Ak);
            isspA(p,k) = (nnzAk < spdensity*n2) || (numblk > 1);
            if ~all(breakyes)
                [I,J] = find(abs(Ak) > 0);
                %%
                %% nonzero elements of Ak.
                %%
                if (breakyes(1) == 0);
                    if (nnzAk <= nztol)
                        idx = find(I<=J);
                        nzlist1(k+1) = nzlist1(k)+length(idx);
                        nzlist2 = [nzlist2; [I(idx), J(idx)] ]; %#ok
                    else
                        nzlist1(k+1:m+1) = inf*ones(1,m-k+1);
                        breakyes(1) = 1;
                    end
                end
                %%
                %% nonzero elements of ||A1||+...+||Ak||.
                %%
                if (breakyes(2) == 0)
                    nztmp = zeros(length(I),1);
                    for t = 1:length(I);
                        i=I(t); j=J(t); nztmp(t)=Asum(i,j);
                    end
                    %% find new nonzero positions when Ak is added to Asum.
                    idx = find(nztmp == 0);
                    nzlist3(k+1) = nzlist3(k) + length(idx);
                    if (nzlist3(k+1) < nztol2);
                        nzlist4 = [ nzlist4; [I(idx), J(idx)] ]; %#ok
                    else
                        nzlist3(k+1:m+1) = inf*ones(1,m-k+1);
                        breakyes(2) = 1;
                    end
                    Asum = Asum+abs(Ak);
                end
            end
        end
        if (numblk == 1)
            isspAy(p,1) = (nzlist1(m+1) < inf) || (nzlist3(m+1) < inf);
        else
            isspAy(p,1) = 1;
        end
        nzlistA{p,1} = nzlist1;    nzlistA{p,2} = nzlist2;
        nzlistAsum{p,1} = nzlist3; nzlistAsum{p,2} = nzlist4;
        %%
        %% nonzero elements of (A1*y1+...Am*ym).
        %%
        if (nzlist3(m+1) < inf);
            if (length(pblk) > 2)
                % m2 = length(pblk{3});
                len = sum(pblk{3});
                DD = spconvert([At{p,3}(:,2:4); len, len, 0]);
                Asum = Asum + abs(At{p,2}*DD*At{p,2}');
            end
            [I,J] = find(Asum > 0);
            if (length(I) < nztol2)
                nzlistAy{p} = [I, J];
            else
                nzlistAy{p} = inf;
            end
        else
            nzlistAy{p} = inf;
        end
    end
end
%%***********************************************************************

author: zarathustr


text:
function [Theta, k] = estimate_plane(X, s)

% [Theta k] = estimate_plane(X)
%
% DESC:
% estimate the parameters of a 3D plane given the pairs [x, y, z]^T
% Theta = [a; b; c; d] where:
%
% a*x1+b*y1+c*z1+d = 0
% a*x2+b*y2+c*z2+d = 0
% a*x3+b*y3+c*z3+d = 0
%
% AUTHOR
% Marco Zuliani - marco.zuliani@gmail.com
%
% VERSION:
% 1.0.0
%
% INPUT:
% X                 = 3D points
% s                 = indices of the points used to estimate the parameter
%                     vector. If empty all the points are used
%
% OUTPUT:
% Theta             = estimated parameter vector Theta = [a; b; c; d]
% k                 = dimension of the minimal subset

% HISTORY:
% 1.0.0             = 07/05/08 - initial version

% cardinality of the MSS
k = 3;

if (nargin == 0) || isempty(X)
    Theta = [];
    return;
end;

if (nargin == 2) && ~isempty(s)
    X = X(:, s);
end;

% check if we have enough points
N = size(X, 2);
if (N < k)
    error('estimate_plane:inputError', ...
        'At least 3 points are required');
end;

A = [transpose(X(1, :)) transpose(X(2, :)) transpose(X(3, :)) ones(N, 1)];
[U S V] = svd(A);
Theta = V(:, 4);

return;

author: RANSAC

repo: RANSAC-Toolbox


text:
function value = biw_condition ( n )

%*****************************************************************************80
%
%% BIW_CONDITION computes the L1 condition of the BIW matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real VALUE, the L1 condition.
%
  if ( n == 1 )
    a_norm = 0.6;
  else
    a_norm = 1.6;
  end

  b_norm = 0.0;
  j = n;
  for i = n : -1 : 1
    aii = 0.5 + i / ( 10 * n );
    if ( i == j )
      bij = 1.0 / aii;
    elseif ( i < j )
      bij = bij / aii;
    end
    b_norm = b_norm + abs ( bij );
  end

  value = a_norm * b_norm;

  return
end

author: johannesgerer


text:
function x = subsref( x, S, cheat )
global cvx___
try
    x = subsref( cvx___.problems( index( x ) ), S );
catch
    error( lasterr );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function class_info=gen_class_info_cityscapes()

class_info=[];

class_info.class_names={'road', 'sidewalk', 'building', 'wall', 'fence', 'pole', 'trafficlight',...
     'trafficsign', 'vegetation', 'terrain', 'sky', 'person', 'rider', 'car', ...
     'truck', 'bus', 'train', 'motorcycle', 'bicycle', 'void'};

class_label_values=uint8([0:18 255]);

class_info.class_label_values=class_label_values;
class_info.background_label_value=uint8(0);
class_info.void_label_values=uint8(255);

cmap=load('cityscape_cmap.mat');
cmap=uint8(cmap.cityscape_cmap);
class_info.mask_cmap=im2double(cmap);

class_info=process_class_info(class_info);

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [Y,XYZmm] = spm_read_vols(V,mask)
% Read in entire image volumes
% FORMAT [Y,XYZmm] = spm_read_vols(V,mask)
% V      - vector of mapped image volumes to read in (from spm_vol)
% mask   - implicit zero mask?
%
% Y      - 4D matrix of image data, fourth dimension indexes images
% XYZmm  - 3xn matrix of XYZ locations returned {mm}
%__________________________________________________________________________
%
% For image data types without a representation of NaN (see spm_type),
% implicit zero masking can be used. If mask is set, then zeros are
% treated as masked, and returned as NaN.
%__________________________________________________________________________
% Copyright (C) 1999-2016 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_read_vols.m 6918 2016-11-02 14:33:11Z guillaume $

%-Argument checks
%--------------------------------------------------------------------------
if nargin<2, mask = false; end

spm_check_orientations(V);

%-Read in image data
%--------------------------------------------------------------------------
n = numel(V);                                           %-#images
Y = zeros([V(1).dim(1:3),n]);                           %-image data matrix

for i=1:n
    for p=1:V(1).dim(3)
        Y(:,:,p,i) = spm_slice_vol(V(i),spm_matrix([0 0 p]),V(i).dim(1:2),0);
    end
end

%-Apply implicit zero mask for image datatypes without a NaNrep
%--------------------------------------------------------------------------
if mask
    %-Work out images without NaNrep
    im = false(n,1);
    for i=1:n, im(i)=~spm_type(V(i).dt(1),'NaNrep'); end
    %-Mask
    Y(Y(:,:,:,im)==0) = NaN;
end

%-Return as 3D matrix if single image
%--------------------------------------------------------------------------
if n==1, Y=Y(:,:,:,1); end

%-Compute XYZmm co-ordinates (if required)
%--------------------------------------------------------------------------
if nargout > 1
    [R,C,P] = ndgrid(1:V(1).dim(1),1:V(1).dim(2),1:V(1).dim(3));
    RCP     = [R(:)';C(:)';P(:)';ones(1,numel(R))];
    XYZmm   = V(1).mat(1:3,:)*RCP;
end


text:
function [ y, m ] = month_carry_english ( y, m )

%*****************************************************************************80
%
%% MONTH_CARRY_ENGLISH carries a year of months on the English calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, the year and month.
%    On output, M is no greater than 12.
%
  while ( 1 )

    months = year_length_months_english ( y );

    if ( m <= months )
      break
    end

    m = m - months;
    y = y + 1;

  end

  return
end

repo: jburkardt-m


text:
function Fx = central_diff( F, x )
% central_diff   Central Difference Gradient
%
% usage:    gradient = central_diff( F, x )
%
% inputs:   F - Values of a scalar function evaluated at x
%               to be differentiated with respect to x
%           x - vector of monotonically increasing coordinates where F is evaluated,
%               or a scalar for evenly spaced coordinates
%
% output:  gradient - numerically evaluated gradient by:
%                     forward difference at the left end;
%                     backward difference at the right end;
%                     central difference in the interior
%
% Written by:    Lt Col Robert A. Canfield
%                Air Force Institute of Technology
%                AFIT/ENY
%                2950 P Street, Bldg. 640
%                WPAFB, OH  45433-7765
%                Robert.Canfield@afit.af.mil
%
% Created:       10/19/00
% Last modified: 10/01/01
%
% Note:     MATLAB's gradient function is not second-order accurate
%           for un-evenly spaced coordinates.
%           This central_diff function uses the correct central difference
%           formula (for interior points) that is second-order accurate
%           for both evenly and un-evenly spaced coordinates.
%           Tested under MATLAB versions 5.2, 5.3.1, and 6.0.
%
% Alternatively, you may patch MATLAB's gradient function by
% replacing the lines...
%   % Take centered differences on interior points
%   if n > 2
%      h = h(3:n) - h(1:n-2);
%      g(2:n-1,:) = (f(3:n,:)-f(1:n-2,:))./h(:,ones(p,1));
%   end
%
% with...
%   % Take centered differences on interior points
%   if n > 2
%      if all(abs(diff(h,2)) < eps) % only use for uniform h (RAC)
%         h = h(3:n) - h(1:n-2);
%         g(2:n-1,:) = (f(3:n,:)-f(1:n-2,:))./h(:,ones(p,1));
%      else   % new formula for un-evenly spaced coordinates (RAC)
%         h = diff(h); h_i=h(1:end-1,ones(p,1)); h_ip1=h(2:end,ones(p,1));
%         g(2:n-1,:) =  (-(h_ip1./h_i).*f(1:n-2,:) + ...
%                         (h_i./h_ip1).*f(3:n,:)   )./ (h_i + h_ip1) + ...
%                         ( 1./h_i - 1./h_ip1 ).*f(2:n-1,:);
%      end
%   end

% Ensure compatible vectors and x monotonically increasing or decreasing
if nargin<2
	m = 1;
   h = 1;
   x = 1;
else
	m = length(x);
end
Tflag = 0;
if ndims(F)==2 & size(F,1)==1    % Treat row vector as a column vector
   F = F.';
   Tflag = 1;
end;

[n,p] = size(F);
if m==0
   error('x cannot be null')
elseif m~=1 & m~=n
	error('First dimension of F and x must be the same.')
elseif m>1 & ~( all(diff(x)>0) | all(diff(x)<0) )
	error('Vector x must be monotonically increasing or decreasing.')
elseif n<=1
   Fx = F / x;
   return
end

Fx = zeros(size(F));
x  = x(:);

% Forward difference at left end
if m==1
	h = x;
else
	h = x(2) - x(1);
end
Fx(1,:) = ( F(2,:) - F(1,:) ) / h;

% Backward difference at right end
if m>1, h = x(m) - x(m-1); end
Fx(n,:) = ( F(n,:) - F(n-1,:) ) / h;

% Central Difference in interior
if n > 3
	if m==1 | all(diff(x)==h) % Evenly spaced formula used in MATLAB's gradient routine
		Fx(2:n-1) = ( F(3:n,:) - F(1:n-2,:) ) / (2*h);
   else
		h = diff(x); h_i=h(1:m-2,ones(p,1)); h_ip1=h(2:m-1,ones(p,1));
		Fx(2:n-1,:) =  (-(h_ip1./h_i).*F(1:n-2,:) + ...
		                 (h_i./h_ip1).*F(3:n,:)   )./ (h_i + h_ip1) + ...
		        ( 1./h_i - 1./h_ip1 ).*F(2:n-1,:);
	end
end
if Tflag, Fx=Fx.'; end

author: Sable

repo: mcbench-benchmarks