text:
function Lrho = construct_connection_Laplacian(measurements)
%function Lrho = construct_connection_Laplacian(measurements)
%
% This function computes and returns the connection Laplacian for the
% rotational observations (see eq. (15) of the paper).

D = length(measurements.t{1}); % D = dimension of SE(d)
N = max(max(measurements.edges));  % N = number of nodes in the pose graph
M = size(measurements.edges,1); % M = number of edges in the pose graph

% The number of nonzero elements in the connection Laplacian; there are
% 2*D^2*M nonzero elements corresponding to the off-diagonal blocks
% containing the raw measurements (one dxd rotation matrix above the main
% diagonal, and another below, for each of the M measurements, plus D*N
% nonzero elements along the main diagonal for the scaled identity
% matrices.  Note that in the case of multiple observations between two
% nodes x_i and x_j the corresponding (i,j) block in the connection
% Laplacian should contain the *sum* of the weighted observations
% kappa_{ij} R_{ij} between these nodes.  We achieve this result below by
% exploiting the fact that MATLAB's sparse() command forms the resulting
% sparse matrix by *summing* values with the same (i,j) indices.

% Copyright (C) 2016 by David M. Rosen

D2 = D^2;
off_diag_inc = 2*D2;

NNZ = off_diag_inc * M + D*N;

%Allocate storage for the row, column, and value vectors

rows = zeros(1, NNZ);
cols = zeros(1, NNZ);
vals = zeros(1, NNZ);

degs = zeros(1, N);  %Vector to store the degrees of the nodes

%Iterate over the measurements in the pose graph
for k = 1:M

    %EXTRACT MEASUREMENT DATA
    i = measurements.edges(k, 1);  %The node that this edge leaves
    j = measurements.edges(k, 2);  %The node that this edge enters

    Rij = measurements.R{k};  %The rotation matrix for this observation

    kappa = measurements.kappa{k};  %The precision for this rotational observation

    %PROCESS MEASUREMENT DATA

    %Increment the degrees of the ith and jth nodes by kappa
    degs(i) = degs(i) + kappa;
    degs(j) = degs(j) + kappa;

    %Set the (i,j)th DxD block of Lrho = -kappa*Rij
    [r, c, Rvect] = rcvize_matrix(Rij, i, j);
    rows(off_diag_inc*(k - 1) + 1 : off_diag_inc*(k-1) + D2) = r;
    cols(off_diag_inc*(k - 1) + 1 : off_diag_inc*(k-1) + D2) = c;
    vals(off_diag_inc*(k - 1) + 1 : off_diag_inc*(k-1) + D2) = -kappa*Rvect;

    %Set the (j,i)th DxD block of Lrho = -kappa*R_ij^T
    [r, c, Rvect] = rcvize_matrix(Rij', j, i);
    rows(off_diag_inc*(k - 1) + D2 + 1 : off_diag_inc*(k-1) + off_diag_inc) = r;
    cols(off_diag_inc*(k - 1) + D2 + 1 : off_diag_inc*(k-1) + off_diag_inc) = c;
    vals(off_diag_inc*(k - 1) + D2 + 1 : off_diag_inc*(k-1) + off_diag_inc) = -kappa*Rvect;
end

%Now set the diagonal elements to be D-fold copies of the degrees of each
%node

rows(off_diag_inc*M + 1 : NNZ) = [1:D*N];
cols(off_diag_inc*M + 1 : NNZ) = [1:D*N];
vals(off_diag_inc*M + 1 : NNZ) = kron(degs, ones(1,D));

Lrho = sparse(rows, cols, vals, D*N, D*N);

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function varargout = smoothMesh(varargin)
%SMOOTHMESH Smooth mesh by replacing each vertex by the average of its neighbors.
%
%   V2 = smoothMesh(V, F)
%   [V2, F2] = smoothMesh(V, F)
%   Performs smoothing of the values given in V, by using adjacency
%   information given in F.
%   V is a numeric array representing either vertex coordinate, or value
%   field associated to each vertex. F is an array of faces, given either
%   as a NF-by-3 or NF-by-4 numeric array, or as a cell array.
%   Artifact adjacencies are added if faces have more than 4 vertices.
%
%   ... = smoothMesh(V, F, NITER)
%   Repeat the smoothing procedure NITER times. This is equivalent to
%   calling the smoothMesh function NITER times.
%
%
%   Example
%     [v f] = torusMesh([50 50 50 30 10 30 45]);
%     v = v + randn(size(v));
%     [v2 f] = smoothMesh(v, f, 3);
%     figure; drawMesh(v2, f);
%     l = light; lighting gouraud
%
%   See also
%     meshes3d, meshAdjacencyMatrix, triangulateFaces, drawMesh
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2013-04-29, using Matlab 7.9.0.529 (R2009b)
% Copyright 2013-2022 INRA - Cepia Software Platform

var1 = varargin{1};
if isstruct(var1)
    vertices = var1.vertices;
    faces = var1.faces;
    varargin(1) = [];
else
    vertices = varargin{1};
    faces = varargin{2};
    varargin(1:2) = [];
end

% determine number of iterations
nIter = 1;
if ~isempty(varargin)
    nIter = varargin{1};
end

% compute adjacency matrix,
% result is a Nv-by-Nv matrix with zeros on the diagonal
adj = meshAdjacencyMatrix(faces);

% ensure the size of the matrix is Nv-by-Nv
% (this can not be the case if some vertices are not referenced)
nv = size(vertices, 1);
if size(adj, 1) < nv
    adj(nv, nv) = 0;
end

% Add "self adjacencies"
adj = adj + speye(nv);

% weight each vertex by the number of its neighbors
w = spdiags(full(sum(adj, 2).^(-1)), 0, nv, nv);
adj = w * adj;

% do averaging to smooth the field
v2 = vertices;
for k = 1:nIter
    v2 = adj * v2;
end

varargout = formatMeshOutput(nargout, v2, faces);

%% Old version
% % Compute vertex adjacencies
% edges = computeMeshEdges(faces);
% v2 = zeros(size(vertices));
%
% % apply several smoothing
% for iter = 1:nIter
%
%     % replace the coords of each vertex by the average coordinate in the
%     % neighborhood
%     for i = 1:size(vertices, 1)
%         edgeInds = sum(edges == i, 2) > 0;
%         neighInds = unique(edges(edgeInds, :));
%         v2(i, :) = mean(vertices(neighInds, :));
%     end
%
%     % update for next iteration
%     vertices = v2;
% end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function dep = cfg_vout_file_fplist(job)

% function dep = cfg_vout_file_fplist(job)
%
% Virtual outputs for cfg_run_file_fplist. Struct with fields .files and
% .dirs. See help on cfg_getfile.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_vout_file_fplist.m 3990 2010-07-13 10:27:17Z volkmar $

rev = '$Rev: 3990 $'; %#ok

dep(1)            = cfg_dep;
if ~ischar(job.filter) || strcmp(job.filter,'<UNDEFINED>')
    dep(1).sname  = 'Selected Files';
else
    dep(1).sname  = sprintf('Selected Files (%s)', job.filter);
end
dep(1).src_output = substruct('.','files');
dep(1).tgt_spec   = cfg_findspec({{'class','cfg_files','strtype','e'}}); % This will display for both file and dir inputs. There is no notion of "no directories" in matching.
dep(2)            = cfg_dep;
dep(2).sname      = 'Subdirectories';
dep(2).src_output = substruct('.','dirs');
dep(2).tgt_spec   = cfg_findspec({{'filter','dir','strtype','e'}});

author: spm


text:
function [params, names] = dexpKernExtractParam(kern)

% DEXPKERNEXTRACTPARAM Extract parameters from the double exponential's
% kernel structure.
%
% FORMAT
% DESC extracts parameters from the double exponential kernel's structure
% into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN params : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC extracts parameters and parameter names from the double exponential
% kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing names for each
% parameter.
%
% SEEALSO dexpKernParamInit, dexpKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : David Luengo, 2009

% KERN

params = [kern.decay kern.variance];
if nargout > 1
  names={'decay', 'variance'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [d dt pred] = bfs_in_mbgl(A,u)
% BFS_IN_MBGL Reimplement the BFS function with MatlabBGL visitors.
%
% [d dt pred] = bfs_in_mbgl(A,u)
%
% See BFS
%
% Example:
%    load ../graphs/bfs_example.mat
%    d = bfs_in_mbgl(A,1)

ip_d = ipdouble(-ones(num_vertices(A),1));
ip_dt = ipdouble(-ones(num_vertices(A),1));
ip_pred = ipdouble(zeros(1,num_vertices(A)));

ip_time = ipdouble(1);

    function time_discover_vertex(u)
        ip_dt(u) = ip_time(1);
        ip_time(1) = ip_time(1) + 1;
    end

    function distance_tree_edge(ei,u,v)
        ip_d(v) = ip_d(u)+1;
    end

    function pred_tree_edge(ei,u,v)
        ip_pred(v) = u;
    end

vis_distance = struct('tree_edge', @distance_tree_edge);
vis_time = struct('discover_vertex', @time_discover_vertex);
vis_pred = struct('tree_edge', @pred_tree_edge);

vis = combine_visitors(vis_distance, vis_time, vis_pred);

ip_d(u) = 0;
ip_dt(u) = ip_time(1);
breadth_first_search(A,u,vis);

d = double(ip_d);
dt = double(ip_dt);
pred = double(ip_pred);

% end the function
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [dat,nvox] = iimg_cluster_extent(dat,volInfo,k)
% Apply a cluster size threshold to a series of index image data vectors
% (each img is one column)
% given volInfo (must be created with iimg_read_img, with extended output)
%
% :Usage:
% ::
%
%     [dat,nvox] = iimg_cluster_extent(dat,volInfo,k)
%
% :Inputs:
%
%   **dat:**
%        may be an indexed image of all image values or only those in mask
%        defined by volInfo

if isempty(k) || k < 2, nvox = []; return, end % do nothing if extent threshold

if size(dat,1) == volInfo.nvox
    dattype = 'full';
elseif size(dat,1) == volInfo.n_inmask
    dattype = 'masked';
else error('data vector does not match size of image in volInfo!')
end

switch dattype
    case 'full', [clindx,nvox] = iimg_cluster_index(dat(volInfo.wh_inmask, :), volInfo.xyzlist', k);
    case 'masked', [clindx,nvox] = iimg_cluster_index(dat,volInfo.xyzlist',k);
end

switch dattype
    case 'full'
        for i = 1:size(dat,2)
            clindx2 = zeros(volInfo.nvox,1);    % put in original index dims
            clindx2(volInfo.wh_inmask) = clindx(:,i);
            dat(:,i) = dat(:,i) .* (clindx2 > 0);
        end

    case 'masked'
        dat = dat .* (clindx > 0);
end

return


text:
classdef MimModelCallback  < CoreBaseClass

    properties (Access = private)
        ModelId
        ModelMap
        Mim
    end

    methods
        function obj = MimModelCallback(modelId, modelMap, mim)
            obj.ModelId = modelId;
            obj.ModelMap = modelMap;
            obj.Mim = mim;
        end

        function value = getModelValue(obj, modelId)
            % Returns the current value of the specified model

           value = obj.ModelMap.getDependentValue(modelId, obj.ModelId);
        end

        function setModelValue(obj, modelId, value)
            % Sets the current value of the specified model

           obj.ModelMap.setValue(modelId, value);
        end

        function addModelItem(obj, modelId, itemId)
           obj.ModelMap.addItem(modelId, itemId);
        end

        function removeModelItem(obj, modelId, itemId)
           obj.ModelMap.removeItem(modelId, itemId);
        end

        function modelId = buildModelId(obj, modelClassName, parameters)
            modelId = obj.ModelMap.buildModelId(modelClassName, parameters);
        end

        function mim = getMim(obj)
            mim = obj.Mim;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%% Radar Plan Position Indicator(PPI) Simulation

% Created by
% R.Vivek
% MIT (Anna University)

%% Description
% This function displays the Radar PPI plot for detection and ...
% tracking of a single target approaching at a specified angle.

%%

function PPI_plot(dist,Angle,Max_Range)

%% "dist" is an array containing the distance (in meters) of the target

%% "Angle" is the angle (in radians) with respect to the radar transceiver
%  at which the target is approaching

%% "Max_Range" is the maximum detectable range of the radar

%% -------------------------------------------------------------------------

angle=repmat(Angle,1,(length(dist)));

for j=1:length(dist)
    for i=(2*pi):-(pi/180):(pi/180)
        [x,y]=(pol2cart(i,55));
        x=ceil(x);
        z=complex(x,y);
        h=compass(z);
        title('Radar PPI Plot');
        ph=findall(gca,'type','patch');
        set(ph,'facecolor',[0,0.5,0]);
        set( h, 'Color', [0.5,1,0] );
        set( h, 'LineWidth', 4 );
        for k = 1:length(h)
            a = get(h(k), 'xdata');
            b = get(h(k), 'ydata');
            set(h(k), 'xdata', a(1:2), 'ydata', b(1:2)')
        end

        if(j>1)
            if(angle(j-1)<i)
                hold on
                plothandle=polar(angle(j-1),dist(j-1),'o');
                set(plothandle,'markerfacecolor',[0,1,0],'markeredgecolor',[0,1,0]);
                get( plothandle ) ;
                set( plothandle, 'Color', [ 0.5, 1, 0 ] );
                set ( plothandle, 'LineWidth', 2 );
                hold off
            end
        end
        if(dist(j)<=Max_Range)
        if(angle(j)>=i)
            hold on
            plothandle=polar(angle(j),dist(j),'o');
            set(plothandle,'markerfacecolor',[0,1,0],'markeredgecolor',[0,1,0])
            get( plothandle ) ;
            set( plothandle, 'Color', [ 0.5, 1, 0 ] );
            set ( plothandle, 'LineWidth', 2 );
            hold off
        end
        end
        pause(1*10^-6);
        end
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function f=maxdiag(x,varargin)
% Maximizaton of the diagonality of hte matrix.
covmat = varargin{1};   % #pix X #pix X #taus
sizevec = varargin{2};
W = reshape(x, sizevec(1)*sizevec(2), sizevec(3)); % #pix X #comp
% W = max(W,eps);
ntau = size(covmat,3); % #taus (different time delays)
ftmp = zeros(1,ntau);
fcol = zeros(1,sizevec(3));
for tau=1:ntau
    covmattmp = covmat(:,:,tau);
    for col=1:sizevec(3)
        wcol = W(:,col);
        fcol(col) = log(wcol'*covmattmp*wcol);
    end
    ftmp(tau) = 0.5*sum(fcol);
end

f = sum(ftmp);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function mv = mrViewVoxelData(ui, roi);
%
% mv = mrViewTimeCourse(<ui>, <roi=selected ROI>);
%
% Plot a multi voxel UI for the selected ROI in a mrViewer UI.
%
% roi can be specified as an index into the mrViewer ROIs, a name
% of the ROI, or a loaded struct. ui defaults to the current viewer.
%
% ras, 07/06.
if notDefined('ui'), ui = mrViewGet; end
if notDefined('roi'),    roi = ui.rois(ui.settings.roi);    end

if ischar(roi)
    roi = ui.rois(cellfind({ui.rois.name}, roi));

elseif isnumeric(roi)
    roi = ui.rois(roi);

end

% this version is for mrVista 2, so we need a loaded mrSESSION.mat file,
% globals and all that (see sessionGUI_viewTimeCourse):
mrGlobals2;

% check if an event-related scan group is assigned.
% If so, use that; otherwise use the GUI selected scans.
scans = GUI.settings.scan;
dt = GUI.settings.dataType;
if isfield(dataTYPES(dt).scanParams(scans(1)), 'scanGroup') & ...
        ~isempty(dataTYPES(dt).scanParams(scans(1)).scanGroup)
    [scans dt] = er_getScanGroup(INPLANE{1});
end

roi = roiParse(roi);
roi.viewType = 'Inplane'; % always in this space, even if it was
                          % loaded from Volume/Gray

mv = multiVoxelUI(INPLANE{1}, roi, scans, dt);

if ~checkfields(GUI, 'controls', 'tcui')
    GUI.controls.tcui = [];
end
GUI.controls.tcui(end+1) = mv.ui.fig;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [funOut, indexStart, problemDom, coeffs, totalDiffOrders] = ...
    toFirstOrder(funIn, rhs, domain, numArgs, cellArg)
%TOFIRSTORDER   Convert higher order anonymous functions to first order systems.
%   Calling sequence:
%      [FUNOUT, INDEXSTART, PROBLEMDOM, COEFFS, TOTALDIFFORDERS] = ...
%          TOFIRSTORDER(FUNIN, RHS, DOMAIN)
%   where the inputs are:
%      FUNIN:  An anonymous function which describes an ODE, usually including
%              higher order derivatives, e.g. @(x,u) diff(u, 2) + sin(u).
%              Usually, the anonymous function comes from the OP field of a
%              CHEBOP.
%      RHS:    The right hand side of the differential equation being solved
%              with a CHEBOP, that is, the right argument of a CHEBOP
%              backslash.
%      DOMAIN: The domain of the problem we're trying to solve.
%   and the outputs are:
%      FUNOUT:     An anonymous function that is the first order reformulation
%                  of FUNIN, which can be passed to the MATLAB solvers.
%      INDEXSTART: A vector that denotes at which index we should start
%                  indexing each variable from so that they're in the correct
%                  order for the MATLAB ODE solvers.
%      PROBLEMDOM: The domain on which the problem can be specified, which may
%                  include breakpoints originally not included in DOMAIN.
%      COEFFS:     A cell array of the coefficients the multiply the highest
%                  order derivatives in problem. For example, for the problem
%                  @(x,u) 5*diff(u, 2) + u, we will have COEFFS{1} = 5.
%      TOTALDIFFORDERS:
%                  A vector that contains the maximum diffOrder applying to
%                  each variable in the problem.
%
%   In addition, one can pass the following arguments to the function, which are
%   only required if working in CHEBMATRIX syntax, for example,
%   @(x,u)[diff(u{1}) + u{2};...]:
%       NVARS:   The number of unknown variables in the ODE.
%       CELLARG: Specify as TRUE if FUNIN is specified in CHEBMATRIX syntax.

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

% Independent variable on the domain
t = chebfun(@(t) t, domain);

% See if we need to determine the number of variables involved, or if it was
% passed as an argument.
if ( nargin < 4)
    % We always have at least one argument to FUNIN, even in the scalar case. In
    % the system case, the first argument to funIn must be the independent time
    % variable. So the number of treeVar arguments needed is one less:
    numArgs = max(1, nargin(funIn) - 1);
    % If working in CHEBMATRIX syntax, we need all possible input variables to
    % be passed in as arguments:
    cellArg = false;
end
args = cell(numArgs, 1);
argsVec = zeros(1, numArgs);

% Populate the args cell
for argCount = 1:numArgs
    argsVec(argCount) = 1;
    args{argCount} = treeVar(argsVec, domain);
    % Reset the index vector
    argsVec = 0*argsVec;
end

% Evaluate FUNIN with the TREEVAR arguments. Need different calling sequences
% depending on whether we're working with CHEBMATRIX syntax and cells or not:
if ( cellArg )
    fevalResult = funIn(t, args);
else
    % If FUNIN only takes one argument, we just give it the TREEVAR argument.
    % Otherwise, the first input will be the independent variable on the domain:
    if ( nargin(funIn) == 1 )
        fevalResult = funIn(args{:});
    else
        fevalResult = funIn(t, args{:});
    end
end

% If we got passed the problem as N\0, i.e. a RHS that did not match the
% dimensions, we need to repmat it so that the RHS variable has the correct
% dimensions below:
if ( isnumeric(rhs) && length(rhs) == 1 )
    rhs = repmat(rhs, size(fevalResult));
end

% Ensure RHS is a CHEBMATRIX
if ( ~isa(rhs, 'chebmatrix') )
    rhs = chebmatrix(rhs, domain);
end

% Initialize cells to store the infix forms of the expressions, the coefficients
% multiplying the highest order derivatives and any variables that appear in the
% anonymous functions:
systemInfix = cell(length(fevalResult), 1);
coeffs = systemInfix;
varArrays = systemInfix;

% We want to return all potential breakpoints caused by piecewise coefficients
% in the problem. So loop through fevalResult, and take the union of the
% breakpoints:
problemDom = fevalResult(1).domain;
for resCounter = 2:length(fevalResult)
    problemDom = union(problemDom, fevalResult(resCounter).domain);
end

% Ensure we also include potential breakpoints from the RHS:
problemDom = union(problemDom, rhs.domain);

% First look at all diffOrders to ensure we start with the correct indices.
% INDEXSTART denotes at which index we should start indexing each variable from.
% E.g., in the coupled system
%   [v'' + w; v + w'']
% we will have v = u(1), v' = u(2), w = u(3), w' = u(4), so INDEXSTART = [1, 3],
% since v and its derivatives starts getting counted at 1, and w and its
% derivatives start getting counted at 3.
%
% The vector INDEXSTARTDER is similar, but here, we also assign the highest
% order derivative of each variable it's own index. This is so that later on, we
% can correctly evaluate the coefficient multiplying the highest order
% derivative in problem. Thus, in the example above, we'd have
%   v = u(1), v' = u(2), v'' = u(3), w = u(4), w' = u(5), w'' = u(6)
% Here, INDEXSTARTDER = [1 4].

% First, we need to find out the maximum diffOrder of each variable appearing in
% the problem. We loop through the each component of the evaluation tree:
totalDiffOrders = zeros(1, numArgs);
for wCounter = 1:length(fevalResult)
    totalDiffOrders = max(totalDiffOrders, ...
        fevalResult(wCounter).tree.diffOrder);
end

% We always start indexing the first variable and its derivative(s) at 1. The
% index of the other variables depend on the cumulative diffOrders of the
% variables with lower indices.
indexStart = [ 1, cumsum(totalDiffOrders(1:end-1)) + 1 ];
% To get the indices of the derivatives as well, we shift all the previous
% indices right by 1, in cumulative fashion:
indexStartDer = indexStart + (0:numArgs-1);

% COEFFARG will be used to evaluate the functions that gives us information
% about the coefficients multiplying the highest order derivative in each
% equation. The vector has to be equally long to the total number of derivatives
% appearing in the problem; we'll then change one of the entries to 1 at a time
% to get the coefficient information.
coeffArg = zeros(1, indexStartDer(end) + totalDiffOrders(end));

% Go through each componenent from the result of evaluating FUNIN,
% and change it to infix format.
for wCounter = 1:length(fevalResult)

    % The current result we're looking at.
    res = fevalResult(wCounter);

    % Current diffOrders
    diffOrders = res.tree.diffOrder;

    % Ensure that we never have the highest derivatives of more than one
    % variable appearing in a single equation:
    if ( sum(totalDiffOrders == diffOrders) > 1 )
        error('CHEBFUN:TREEVAR:toFirstOrder:diffOrders', ...
            ['The highest order derivative of more than one variable ' ...
            'appears to be\npresent in the same equation. ' ...
            'Unable to convert to first order format.'])
    end

    % Expand the tree, so that PLUS rather than TIMES is sitting at the top of
    % it.
    expTree = treeVar.expandTree(res.tree, totalDiffOrders);

    % Split the tree into derivative part and non-derivative part.
    [newTree, derTree] = treeVar.splitTree(expTree, totalDiffOrders);

    % If newTree is empty, we only have a derivative part in the expression,
    % e.g. diff(u) = 0. We must replace it with a 0, as otherwise, we can't
    % evaluate the resulting odeFun in the ODE solvers.
    if ( isempty(newTree) )
        newTree = 0;
    end

    % Find what argument corresponds to the highest derivative one in the
    % current expression we're looking at. This will also be the order in which
    % we store the outputs from converting individual expressions to first order
    % form -- if the input is of the form @(x,u,v) [diff(v) - u; diff(u) - v],
    % the equations have to be sorted so that they'll be correctly converted.
    maxDerLoc = find(expTree.diffOrder == totalDiffOrders);

    % Convert the derivative part to infix form.
    % Indicate that we are converting a coeffFun
    isCoeffFun = true;
    [infixDer, varArrayDer] = ...
        treeVar.tree2infix(derTree, maxDerLoc, indexStartDer, isCoeffFun);

    % Convert the infix form of the expression that gives us the coefficient
    % multiplying the highest order derivative appearing in the expression to an
    % anonymous function we can evaluate:
    coeffFun = treeVar.toAnon(infixDer, varArrayDer);

    % Reset coeffArg for next evaluation:
    coeffArg = 0*coeffArg;

    % Replace one of the 0s in coeffFun with 1 so that we can evaluate COEFFFUN:
    if ( maxDerLoc == numArgs )
        % The last variable in the system currently appears in the highest order
        % derivate.
        coeffArg(end) = 1;
    else
        % The variable with index maxDerLoc+1 is the next variable we need to
        % start numbering at. So subtract 1 for the index of the highest
        % derivative we're currently interested in.
        coeffArg(indexStartDer(maxDerLoc+1) - 1) = 1;
    end

    % Evaluate the COEFFFUN to the coefficient!
    coeffs{maxDerLoc} = coeffFun(t, coeffArg);

    % Now work with the remaining syntax tree of the current expression of
    % interest. We need to negate the syntax tree as we're moving it to the
    % right-hand side. But if it already starts with a unary minus, we can
    % simply remove it rather than doing a double negation:
    % [TODO: Remove double UMINUS]
    newTree = struct('method', 'minus', 'numArgs', 2, ...
        'left', rhs{wCounter}, 'right', newTree);
    % Convert current expression to infix form:
    isCoeffFun = false;
    [infix, varArray] = ...
        treeVar.tree2infix(newTree, maxDerLoc, indexStart, isCoeffFun);
    % Store the infix form and the variables that appeared in the anonymous
    % function.
    systemInfix{maxDerLoc} = infix;
    varArrays{maxDerLoc} = varArray;
end

% Convert all the infix expressions, coefficients and variables stored to an
% anonymous function we can evaluate and use as the RHS of our ODE system:
funOut = treeVar.toRHS(systemInfix, varArrays, coeffs, ...
    indexStart, totalDiffOrders);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function a = r8sto_indicator ( n )

%*****************************************************************************80
%
%% R8STO_INDICATOR sets up a R8STO indicator matrix.
%
%  Discussion:
%
%    The R8STO storage format is used for a symmetric Toeplitz matrix.
%    It stores the N elements of the first row.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Output, real A(N), the R8STO matrix.
%
  fac = 10^( i4_log_10 ( n ) + 1 );

  i = 1;
  k = 0;
  for j = 1 : n
    k = k + 1;
    a(k) = fac * i + j;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% test for OpengJPEG implementation of the JPEG2000 Codec
clear options;
name = 'lena';
M = load_image(name);

test_type = 'rate';
test_type = 'db';
test_type = 'bits';

switch test_type
    case 'rate'
        % compression of a factor 30
        options.rate = 30;
    case 'db'
        % compression with given erro
        options.db = 30;
    case 'bits'
        options.nbr_bits = 5e4;
end

% Compression:
[stream,nbr_bits] = perform_jp2k_compression(M,options);
% De-compression:
M1 = perform_jp2k_compression(stream,options);
bd = ceil(log2(max(M(:))));
nbr_bits_orig = bd*prod(size(M));

clf;
subplot(1,2,1);
imagesc(M); axis image, axis off;
title('Original');
subplot(1,2,2);
imagesc(M1); axis image, axis off;
title(['Compressed, db=' num2str(psnr(M,M1,2^bd),3) ', rate=' num2str(nbr_bits_orig/nbr_bits,3), ' nbrbits=' num2str(nbr_bits,3) ]);
colormap gray(256)

author: gpeyre

repo: matlab-toolboxes


text:
function T = slfld(X, nums, varargin)
%SLFLD Performs Fisher Linear Discriminant Analysis
%
% $ Syntax $
%   - T = slfld(X, nums)
%   - T = slfld(X, nums, ...)
%
% $ Arguments $
%   - X:        the training sample matrix
%   - nums:     the numbers of samples in all classes
%   - T:        the solved transform matrix
%
% $ Description $
%   - T = slfld(X, nums) performs Fisher linear discriminant analysis
%     on the samples X in default way.
%
%   - T = slfld(X, nums, ...) performs Fisher linear discriminant analysis
%     on the samples X according to the specified properties.
%     \*
%     \t   Table 1.  The properties of Fisher Discriminant Analysis   \\
%     \h     name     &     description                                \\
%           'prepca'  &  Whether to perform a preamble PCA to first
%                        reduce the dimensions to the samples' rank.
%                        default = false.                              \\
%           'whiten' &  The cell containing the arguments for computing
%                       the whitening transform.  default = {}.       \\
%           'dimset'  &  The cell containing the arguments for determining
%                        the output feature dimension. default = {}.
%                        (refer to sldim_by_eigval).                   \\
%           'Sb'      &  The pre-computed between-class scattering matrix
%                        or the cell containing the arguments for
%                        computing the scatter matrix in the form
%                        {type, ...}, which is input to slscatter.     \\
%           'Sw'      &  The pre-computed within-class scattering matrix
%                        or the cell containing the arguments for
%                        computing the scatter matrix in the form
%                        {type, ...}, which is input to slscatter.     \\
%         'weights'   &  The sample weights. default = [].             \\
%     \*
%
% $ Remarks $
%   -# The function solves the transform in mainly two stages:
%      First, whiten the samples so that the between-class scattering
%      become the identity matrix; then solves the projection to maximize
%      the whitened between-class scattering. If prepca is set to true
%      a preamble step for reducing the dimensions to rank by PCA is taken.
%
%   -# The function is very flexible. By tuning the arguments for
%      pre-pca, whitening, and the computation of Sb and Sw, it turns to
%      be various LDA-based algorithms, including fisher LDA,
%      regularized LDA, weighted pairwise LDA, dual-space LDA etc.
%
%   -# If Sw or its computing rule is given, the whiten transform is
%      directly solved from Sw,  otherwise the whitening transform is
%      solved from samples. If Sb is given, the whitened between-class
%      scattering is computed by directly applying the whiten transform
%      to Sb, otherwise, Sb is computed from whitenned samples. If both
%      Sb and Sw are given, then the samples are not used in the function.
%      In this cases, you can simply input an empty X.
%
%   -# If both Sb and Sw are given, the pre-pca step will not be conducted.
%      no matter whether prepca is true or false.
%
% $ History $
%   - Created by Dahua Lin on Apr 30th, 2006
%   - Modified by Dahua Lin on Sep 10th, 2006
%       - replace sladd by sladdvec to increase efficiency.
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('slfld', 2);
end

% check size

if ~isempty(X)
    if ndims(X) ~= 2
        error('sltoolbox:invaliddims', ...
            'The sample matrix X should be a 2D matrix');
    end
    [d, n] = size(X);

    k = length(nums);
    if ~isequal(size(nums), [1, k]);
        error('sltoolbox:invaliddims', ...
            'The nums vector should be a row vector');
    end
    if sum(nums) ~= n
        error('sltoolbox:sizmismatch', ...
            'The total number in nums is not consistent with that in X');
    end
end

% check options

opts.prepca = false;
opts.whiten = {};
opts.dimset = {};
opts.Sb = {'Sb'};
opts.Sw = {'Sw'};
opts.weights = [];
opts = slparseprops(opts, varargin{:});

has_Sb = ~isempty(opts.Sb) && isnumeric(opts.Sb);
has_Sw = ~isempty(opts.Sw) && isnumeric(opts.Sw);
if has_Sb && has_Sw
    use_samples = false;
    d = size(opts.Sw, 1);

    if ~isequal(size(opts.Sb), [d, d]) || ~isequal(size(opts.Sw), [d, d])
        error('sltoolbox:sizmismatch', ...
            'Size consistency in Sb and Sw');
    end

else
    if isempty(X)
        error('sltoolbox:invalidargs', ...
            'The samples cannot be empty when Sb or Sw is not pre-computed');
    end
    use_samples = true;
    if (has_Sb && ~isequal(size(opts.Sb), [d, d])) || (has_Sw && ~isequal(size(opts.Sw), [d, d]))
        error('sltoolbox:sizmismatch', ...
            'Size consistency in Sb and Sw');
    end

end
w = opts.weights;

%% Compute

%% Step 0: Pre-PCA
pca_computed = false;
if use_samples && opts.prepca
    SPCA = slpca(X, 'weights', w);
    X = SPCA.P' * sladdvec(X, -SPCA.vmean, 1);
    pca_computed = true;
end

%% Step 1: Compute whiten transform

if has_Sw
    TW = slwhiten_from_cov(opts.Sw, opts.whiten{:});
elseif ~isempty(opts.Sw) && ~isequal(opts.Sw, {'Sw'})
    Sw = slscatter(X, opts.Sw{:}, 'sweights', w, 'nums', nums);
    TW = slwhiten_from_cov(Sw);
    clear Sw;
else
    TW = slwhiten_from_samples(make_withinclass_diffvecs(X, w, nums), ...
        'weights', w, opts.whiten{:});
end

if pca_computed
    T1 = SPCA.P * TW;
    clear SPCA TW;
else
    T1 = TW;
    clear TW;
end

%% Step 2: Compute the second-stage transform

if has_Sb
    WSb = T1' * opts.Sb * T1;
else
    X = T1' * X;
    WSb = slscatter(X, opts.Sb{:}, 'sweights', w, 'nums', nums);
end
[evs, T2] = slsymeig(WSb);
rk2 = sldim_by_eigval(evs, opts.dimset{:});
T2 = T2(:, 1:rk2);

%% Integrate the transforms

T = T1 * T2;

%% The function for making the difference vectors
function Y = make_withinclass_diffvecs(X, w, nums)

mvs = slmeans(X, w, nums);
Y = X;
[sp, ep] = slnums2bounds(nums);
k = length(nums);
for i = 1 : k
    Y(:, sp(i):ep(i)) = sladdvec(X(:, sp(i):ep(i)), -mvs(:,i), 1);
end

author: lmthang


text:
function obj = updateSurfc(obj, dataIndex)

    if strcmpi(obj.State.Plot(dataIndex).Class, 'surface')
        surfaceIndex = dataIndex;
        updateSurfOnly(obj, surfaceIndex)
    elseif strcmpi(obj.State.Plot(dataIndex).Class, 'contour')
        contourIndex = dataIndex;
        updateContourOnly(obj, contourIndex)
    end

end

function updateContourOnly(obj, contourIndex)

    %-AXIS INDEX-%
    axIndex = obj.getAxisIndex(obj.State.Plot(contourIndex).AssociatedAxis);

    %-CHECK FOR MULTIPLE AXES-%
    [xsource, ysource] = findSourceAxis(obj,axIndex);

    %-AXIS DATA STRUCTURE-%
    axisData = get(obj.State.Plot(contourIndex).AssociatedAxis);

    %-CONTOUR DATA STRUCTURE- %
    contourData = get(obj.State.Plot(contourIndex).Handle);
    surfData = get(obj.State.Plot(contourIndex-1).Handle);
    figureData = get(obj.State.Figure.Handle);

    %-------------------------------------------------------------------------%

    %-associate scene-%
    obj.data{contourIndex}.scene = sprintf('scene%d', xsource);

    %-------------------------------------------------------------------------%

    %-scatter3d type for contour projection-%
    obj.data{contourIndex}.type = 'scatter3d';
    obj.data{contourIndex}.mode = 'lines';

    %-------------------------------------------------------------------------%

    %-get colormap-%
    cMap = figureData.Colormap;
    fac = 1/(length(cMap)-1);
    colorScale = {};

    for c = 1: length(cMap)
        colorScale{c} = { (c-1)*fac , sprintf('rgb(%f,%f,%f)', 255*cMap(c, :))};
    end

    %-------------------------------------------------------------------------%

    %-get plot data-%
    contourMatrix = contourData.ContourMatrix;

    xData = [];
    yData = [];
    zData = [];
    cData = [];

    zmin = axisData.ZLim(1);
    len = size(contourMatrix, 2);
    n = 1;

    while (n < len)

        %---------------------------------------------------------------------%

        %-get plot data-%
        m = contourMatrix(2, n);
        zlevel = contourMatrix(1, n);

        xData = [xData, contourMatrix(1, n+1:n+m), NaN];
        yData = [yData, contourMatrix(2, n+1:n+m), NaN];
        zData = [zData, zmin * ones(1, m), NaN];

        %---------------------------------------------------------------------%

        %-get edge color-%
        if isnumeric(contourData.LineColor)
            cData = sprintf('rgb(%f,%f,%f)', 255*contourData.LineColor);

        elseif strcmpi(contourData.LineColor, 'interp')
            cData = zData;
            obj.data{contourIndex}.line.colorscale = colorScale;

        elseif strcmpi(contourData.LineColor, 'flat')
            [err, r] = min(abs(surfData.ZData - zlevel));
            [~, c] = min(err);
            r = r(c);

            cData = [cData, surfData.ZData(r, c) * ones(1, m), NaN];
            obj.data{contourIndex}.line.colorscale = colorScale;

        elseif strcmpi(contourData.LineColor, 'none')
            cData = 'rgba(0,0,0,0)';

        end

        n = n + m + 1;

    end

    %-------------------------------------------------------------------------%

    %-set data on scatter3d-%
    obj.data{contourIndex}.x = xData;
    obj.data{contourIndex}.y = yData;
    obj.data{contourIndex}.z = zData;

    %-------------------------------------------------------------------------%

    %-set edge color-%
    obj.data{contourIndex}.line.color = cData;

    %-------------------------------------------------------------------------%

    %-line style-%

    obj.data{contourIndex}.line.width = 2*contourData.LineWidth;

    switch contourData.LineStyle
        case '-'
            obj.data{contourIndex}.line.dash = 'solid';
        case '--'
            obj.data{contourIndex}.line.dash = 'dash';
        case '-.'
            obj.data{contourIndex}.line.dash = 'dashdot';
        case ':'
            obj.data{contourIndex}.line.dash = 'dot';
    end

    %-------------------------------------------------------------------------%

    %-surface name-%
    obj.data{contourIndex}.name = contourData.DisplayName;

    %-------------------------------------------------------------------------%

    %-surface showscale-%
    obj.data{contourIndex}.showscale = false;

    %-------------------------------------------------------------------------%

    %-surface visible-%
    obj.data{contourIndex}.visible = strcmp(contourData.Visible,'on');

    %-------------------------------------------------------------------------%
end

function updateSurfOnly(obj, surfaceIndex)

    %-AXIS INDEX-%
    axIndex = obj.getAxisIndex(obj.State.Plot(surfaceIndex).AssociatedAxis);

    %-CHECK FOR MULTIPLE AXES-%
    [xsource, ysource] = findSourceAxis(obj,axIndex);

    %-SURFACE DATA STRUCTURE- %
    meshData = get(obj.State.Plot(surfaceIndex).Handle);
    figureData = get(obj.State.Figure.Handle);

    %-AXIS STRUCTURE-%
    axisData = get(ancestor(meshData.Parent,'axes'));

    %-SCENE DATA-%
    eval( sprintf('scene = obj.layout.scene%d;', xsource) );

    %-GET CONTOUR INDEX-%
    obj.PlotOptions.nPlots = obj.PlotOptions.nPlots + 1;
    contourIndex = obj.PlotOptions.nPlots;
    obj.PlotOptions.contourIndex(surfaceIndex) = contourIndex;

    %-------------------------------------------------------------------------%

    %-associate scene-%
    obj.data{surfaceIndex}.scene = sprintf('scene%d', xsource);
    obj.data{contourIndex}.scene = sprintf('scene%d', xsource);

    %-------------------------------------------------------------------------%

    %-surface type for face color-%
    obj.data{surfaceIndex}.type = 'surface';

    %-scatter3d type for contour mesh lines-%
    obj.data{contourIndex}.type = 'scatter3d';
    obj.data{contourIndex}.mode = 'lines';

    %-------------------------------------------------------------------------%

    %-get plot data-%
    xData = meshData.XData;
    yData = meshData.YData;
    zData = meshData.ZData;

    if isvector(xData)
        [xData, yData] = meshgrid(xData, yData);
    end

    %-reformat data to mesh-%
    xDataSurface = xData;
    yDataSurface = yData;
    zDataSurface = zData;

    xDataContourDir1 = [xDataSurface; NaN(1, size(xDataSurface, 2))];
    yDataContourDir1 = [yDataSurface; NaN(1, size(yDataSurface, 2))];
    zDataContourDir1 = [zDataSurface; NaN(1, size(zDataSurface, 2))];

    xDataContourDir2 = xDataContourDir1(1:end-1,:)';
    yDataContourDir2 = yDataContourDir1(1:end-1,:)';
    zDataContourDir2 = zDataContourDir1(1:end-1,:)';

    xDataContourDir2 = [xDataContourDir2; NaN(1, size(xDataContourDir2, 2))];
    yDataContourDir2 = [yDataContourDir2; NaN(1, size(yDataContourDir2, 2))];
    zDataContourDir2 = [zDataContourDir2; NaN(1, size(zDataContourDir2, 2))];

    xDataContour = [xDataContourDir1(:); xDataContourDir2(:)];
    yDataContour = [yDataContourDir1(:); yDataContourDir2(:)];
    zDataContour = [zDataContourDir1(:); zDataContourDir2(:)];

    %-------------------------------------------------------------------------%

    %-set data on surface-%
    obj.data{surfaceIndex}.x = xDataSurface;
    obj.data{surfaceIndex}.y = yDataSurface;
    obj.data{surfaceIndex}.z = zDataSurface;

    %- setting grid mesh by default -%
    % x-direction
    xData = xData(1, :);
    obj.data{surfaceIndex}.contours.x.start = xData(1);
    obj.data{surfaceIndex}.contours.x.end = xData(end);
    obj.data{surfaceIndex}.contours.x.size = mean(diff(xData));
    obj.data{surfaceIndex}.contours.x.show = true;

    % y-direction
    yData = yData(:, 1);
    obj.data{surfaceIndex}.contours.y.start = yData(1);
    obj.data{surfaceIndex}.contours.y.end = yData(end);
    obj.data{surfaceIndex}.contours.y.size = mean(diff(yData));;
    obj.data{surfaceIndex}.contours.y.show = true;

    %-------------------------------------------------------------------------%

    %-set data on scatter3d-%
    obj.data{contourIndex}.x = xDataContour(:);
    obj.data{contourIndex}.y = yDataContour(:);
    obj.data{contourIndex}.z = zDataContour(:);

    %-------------------------------------------------------------------------%

    %-COLORING-%

    %-------------------------------------------------------------------------%

    %-get colormap-%
    cMap = figureData.Colormap;
    fac = 1/(length(cMap)-1);
    colorScale = {};

    for c = 1: length(cMap)
        colorScale{c} = { (c-1)*fac , sprintf('rgb(%f,%f,%f)', 255*cMap(c, :))};
    end

    %-------------------------------------------------------------------------%

    %-get edge color-%
    if isnumeric(meshData.EdgeColor)
        cDataContour = sprintf('rgb(%f,%f,%f)', 255*meshData.EdgeColor);

    elseif strcmpi(meshData.EdgeColor, 'interp')
        cDataContour = zDataContour(:);
        obj.data{contourIndex}.line.colorscale = colorScale;
        obj.data{surfaceIndex}.contours.x.colorscale = cDataContour;
        obj.data{surfaceIndex}.contours.y.colorscale = cDataContour;

        obj.data{surfaceIndex}.contours.x.show = false;
        obj.data{surfaceIndex}.contours.y.show = false;

    elseif strcmpi(meshData.EdgeColor, 'flat')
        cData = meshData.CData;

        if size(cData, 3) ~= 1
            cMap = unique( reshape(cData, ...
                [size(cData,1)*size(cData,2), size(cData,3)]), 'rows' );
            cData = rgb2ind(cData, cMap);

            edgeColorScale = {};
            fac = 1/(length(cMap)-1);

            for c = 1: length(cMap)
                edgeColorScale{c} = { (c-1)*fac , sprintf('rgb(%f,%f,%f)', 255*cMap(c, :))};
            end

            obj.data{surfaceIndex}.line.cmin = 0;
            obj.data{surfaceIndex}.line.cmax = 255;
            obj.data{contourIndex}.line.colorscale = edgeColorScale;
        else
            obj.data{contourIndex}.line.cmin = axisData.CLim(1);
            obj.data{contourIndex}.line.cmax = axisData.CLim(2);
            obj.data{contourIndex}.line.colorscale = colorScale;
        end

        cDataContourDir1 = [cData; NaN(1, size(cData, 2))];
        cDataContourDir2 = cDataContourDir1(1:end-1,:)';
        cDataContourDir2 = [cDataContourDir2; NaN(1, size(cDataContourDir2, 2))];
        cDataContour = [cDataContourDir1(:); cDataContourDir2(:)];

        obj.data{surfaceIndex}.contours.x.show = false;
        obj.data{surfaceIndex}.contours.y.show = false;

    elseif strcmpi(meshData.EdgeColor, 'none')
        cDataContour = 'rgba(0,0,0,0)';

    end

    %-set edge color-%
    obj.data{contourIndex}.line.color = cDataContour;
    obj.data{surfaceIndex}.contours.x.color = cDataContour;
    obj.data{surfaceIndex}.contours.y.color = cDataContour;

    %-------------------------------------------------------------------------%

    %-get face color-%
    faceColor = meshData.FaceColor;

    if isnumeric(faceColor)

        if all(faceColor == [1, 1, 1])
            faceColor = [0.96, 0.96, 0.96];
        end

        for n = 1:size(zDataSurface, 2)
            for m = 1:size(zDataSurface, 1)
                cDataSurface(m, n, :) = faceColor;
            end
        end

        [cDataSurface, cMapSurface] = rgb2ind(cDataSurface, 256);
        cDataSurface = double(cDataSurface) + axisData.CLim(1);

        for c = 1: size(cMapSurface, 1)
            colorScale{c} = { (c-1)*fac , sprintf('rgba(%f,%f,%f, 1)', cMapSurface(c, :))};
        end

        obj.data{surfaceIndex}.cmin = axisData.CLim(1);
        obj.data{surfaceIndex}.cmax = axisData.CLim(2);

    elseif strcmpi(faceColor, 'interp')
        cDataSurface = zDataSurface;

        if surfaceIndex > xsource
            cData = [];

            for idx = xsource:surfaceIndex
                cData = [cData; obj.data{idx}.z];
            end

            cMin = min(cData(:));
            cMax = max(cData(:));

            for idx = xsource:surfaceIndex
                obj.data{idx}.cmin = cMin;
                obj.data{idx}.cmax = cMax;
            end
        end

    elseif strcmpi(faceColor, 'flat')
        cData = meshData.CData;

        if size(cData, 3) ~= 1
            cMap = unique( reshape(cData, ...
                [size(cData,1)*size(cData,2), size(cData,3)]), 'rows' );
            cDataSurface = rgb2ind(cData, cMap);

            colorScale = {};
            fac = 1/(length(cMap)-1);

            for c = 1: length(cMap)
                colorScale{c} = { (c-1)*fac , sprintf('rgb(%f,%f,%f)', 255*cMap(c, :))};
            end
        else
            cDataSurface = cData;
        end

    end

    %-set face color-%
    obj.data{surfaceIndex}.colorscale = colorScale;
    obj.data{surfaceIndex}.surfacecolor = cDataSurface;

    %-lighting settings-%

    if isnumeric(meshData.FaceColor) && all(meshData.FaceColor == [1, 1, 1])
        obj.data{surfaceIndex}.lighting.diffuse = 0.5;
        obj.data{surfaceIndex}.lighting.ambient = 0.725;
    else
        % obj.data{surfaceIndex}.lighting.diffuse = 1.0;
        % obj.data{surfaceIndex}.lighting.ambient = 0.9;
    end

    if meshData.FaceAlpha ~= 1
        obj.data{surfaceIndex}.lighting.diffuse = 0.5;
        obj.data{surfaceIndex}.lighting.ambient = 0.725 + (1-meshData.FaceAlpha);
    end

    if obj.PlotlyDefaults.IsLight
        obj.data{surfaceIndex}.lighting.diffuse = 1.0;
        obj.data{surfaceIndex}.lighting.ambient = 0.3;
    end

    %-opacity-%
    obj.data{surfaceIndex}.opacity = meshData.FaceAlpha;

    %-------------------------------------------------------------------------%

    %-line style-%
    obj.data{contourIndex}.line.width = 3*meshData.LineWidth;

    if strcmpi(meshData.LineStyle, '-')
        obj.data{contourIndex}.line.dash = 'solid';
    else
        obj.data{contourIndex}.line.dash = 'dot';
        obj.data{surfaceIndex}.contours.x.show = false;
        obj.data{surfaceIndex}.contours.y.show = false;
    end

    %-------------------------------------------------------------------------%

    %-SCENE CONFIGUTATION-%

    %-------------------------------------------------------------------------%

    %-aspect ratio-%
    asr = obj.PlotOptions.AspectRatio;

    if ~isempty(asr)
        if ischar(asr)
            scene.aspectmode = asr;
        elseif isvector(ar) && length(asr) == 3
            xar = asr(1);
            yar = asr(2);
            zar = asr(3);
        end
    else

        %-define as default-%
        xar = max(xData(:));
        yar = max(yData(:));
        xyar = max([xar, yar]);
        zar = 0.7*xyar;
    end

    scene.aspectratio.x = 1.15*xyar;
    scene.aspectratio.y = 1.0*xyar;
    scene.aspectratio.z = zar;

    %---------------------------------------------------------------------%

    %-camera eye-%
    ey = obj.PlotOptions.CameraEye;

    if ~isempty(ey)
        if isvector(ey) && length(ey) == 3
            scene.camera.eye.x = ey(1);
            scene.camera.eye.y = ey(2);
            scene.camera.eye.z = ey(3);
        end
    else

        %-define as default-%
        xey = - xyar; if xey>0 xfac = 0.1; else xfac = -0.1; end
        yey = - xyar; if yey>0 yfac = -0.5; else yfac = 0.5; end
        if zar>0 zfac = 0.1; else zfac = -0.1; end

        scene.camera.eye.x = xey + xfac*xey;
        scene.camera.eye.y = yey + yfac*yey;
        scene.camera.eye.z = zar + zfac*zar;
    end

    %-------------------------------------------------------------------------%

    %-scene axis configuration-%

    scene.xaxis.range = axisData.XLim;
    scene.yaxis.range = axisData.YLim;
    scene.zaxis.range = axisData.ZLim;

    scene.xaxis.tickvals = axisData.XTick;
    scene.xaxis.ticktext = axisData.XTickLabel;

    scene.yaxis.tickvals = axisData.YTick;
    scene.yaxis.ticktext = axisData.YTickLabel;

    scene.zaxis.tickvals = axisData.ZTick;
    scene.zaxis.ticktext = axisData.ZTickLabel;

    scene.xaxis.zeroline = false;
    scene.yaxis.zeroline = false;
    scene.zaxis.zeroline = false;

    scene.xaxis.showline = true;
    scene.yaxis.showline = true;
    scene.zaxis.showline = true;

    scene.xaxis.tickcolor = 'rgba(0,0,0,1)';
    scene.yaxis.tickcolor = 'rgba(0,0,0,1)';
    scene.zaxis.tickcolor = 'rgba(0,0,0,1)';

    scene.xaxis.ticklabelposition = 'outside';
    scene.yaxis.ticklabelposition = 'outside';
    scene.zaxis.ticklabelposition = 'outside';

    scene.xaxis.title = axisData.XLabel.String;
    scene.yaxis.title = axisData.YLabel.String;
    scene.zaxis.title = axisData.ZLabel.String;

    scene.xaxis.tickfont.size = axisData.FontSize;
    scene.yaxis.tickfont.size = axisData.FontSize;
    scene.zaxis.tickfont.size = axisData.FontSize;

    scene.xaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);
    scene.yaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);
    scene.zaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);

    %-------------------------------------------------------------------------%

    %-SET SCENE TO LAYOUT-%
    obj.layout = setfield(obj.layout, sprintf('scene%d', xsource), scene);

    %-------------------------------------------------------------------------%

    %-surface name-%
    obj.data{surfaceIndex}.name = meshData.DisplayName;
    obj.data{contourIndex}.name = meshData.DisplayName;

    %-------------------------------------------------------------------------%

    %-surface showscale-%
    obj.data{surfaceIndex}.showscale = false;
    obj.data{contourIndex}.showscale = false;

    %-------------------------------------------------------------------------%

    %-surface visible-%
    obj.data{surfaceIndex}.visible = strcmp(meshData.Visible,'on');
    obj.data{contourIndex}.visible = strcmp(meshData.Visible,'on');

    %-------------------------------------------------------------------------%

    leg = get(meshData.Annotation);
    legInfo = get(leg.LegendInformation);

    switch legInfo.IconDisplayStyle
        case 'on'
            showleg = true;
        case 'off'
            showleg = false;
    end

    obj.data{surfaceIndex}.showlegend = showleg;

    %-------------------------------------------------------------------------%
end

author: plotly


text:
function test_int_test04 ( )

%*****************************************************************************80
%
%% TEST04 applies a composite Gauss-Legendre rule to finite interval 1D problems.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  Use a composite 4 point Gauss-Legendre rule,\n' );
  fprintf ( 1, '  for 1D finite interval problems.\n' );

  prob_num = p00_prob_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Problem       Exact\n' );
  fprintf ( 1, '         Ints   Approx	   Error\n' );
%
%  Pick a problem.
%
  for prob = 1 : prob_num

    exact = p00_exact ( prob );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %4d        %14f\n', prob, exact );
%
%  Pick a number of subintervals.
%
    for int_log = 0 : 7

      int_num = 2^int_log;

      result = p00_gauss_legendre ( prob, int_num );

      error = abs ( exact - result );

      fprintf ( 1, '        %4d  %14f  %14e\n', int_num, result, error );

    end

  end

  return
end


text:
%% select file
if ~exist('zip_file_path', 'var') || ~exist(zip_file_path, 'file')
    [file_name, folder_path] = uigetfile('*.zip');
    zip_file_path = fullfile(folder_path, file_name);
end

%% unzip the zip file
[s_dir, s_name, ~] = fileparts(zip_file_path);
output_folder = fullfile(s_dir, s_name);
unzip(zip_file_path, output_folder);

%% select the mat file
temp = dir(output_folder);
for m=1:length(temp)
   [~, ~, file_type] = fileparts(temp(m).name);
   if strcmpi(file_type, '.mat')
      matfile_name = fullfile(output_folder, temp(m).name);
   elseif strcmpi(file_type, '.txt')
       textfile_name = fullfile(output_folder, temp(m).name);
   end
end
if ~exist(matfile_name, 'file')
    error('no mat file in the selected zip file');
end

%% load mat file
load(matfile_name);

% check whether we used the same computer again
if exist(neuron.P.log_file, 'file')
    movefile(matfile_name, neuron.P.log_folder);
    movefile(textfile_name, neuron.P.log_file);
    rmdir(output_folder);
else
    original_log_folder = neuron.P.log_folder;
    original_logfile = neuron.P.log_file;
    neuron.P.log_folder = output_folder;
    neuron.P.log_file = textfile_name;
end

fprintf('The results have been loaded. \n');

author: zhoupc


text:
function [mCatalog] = import_start(sFilterDir,FilePath)
    %GUI file handler (used for MapSeis)
    %13.04.14: small fix allowing to pass in the filepath directly to avoid

    if nargin<2
        FilePath=[];
    end

    % Default return value
    mCatalog = nan;

    % Check if filter directory ends with the appropriate directory separator
    sLastChar = sFilterDir(length(sFilterDir));
    if ~(sLastChar == filesep)
        sFilterDir = [sFilterDir filesep];
    end

    % Get directory in filter directory
    vrDir = dir(sFilterDir);
    if isempty(vrDir)
        errordlg('Directory does not exist');
        return;
    elseif length(vrDir) == 2
        errordlg('Directory is empty. No Filters found.');
        return;
    end

    % Initialize filter list
    cFilter = {};
    % Counter for detected filters
    nFilter = 1;

    for nCnt = 1:length(vrDir)
        % Check for directories
        if ~((strcmp(vrDir(nCnt).name, '.')) || (strcmp(vrDir(nCnt).name, '..')))
            % Get name of file
            sName = vrDir(nCnt).name;
            nPos = strfind(sName, '.m');
            % Is it a .m-file
            if ~isempty(nPos)
                sName = sName(1:(nPos - 1));
                try
                    % Is it an import filter
                    sDescription = feval(sName, 0);
                catch
                    sDescription = [];
                end
                % Is it really an import filter
                if ischar(sDescription)
                    % Insert it into filter list
                    cFilter{nFilter, 1} = sDescription;
                    cFilter{nFilter, 2} = sName;
                    nFilter = nFilter + 1;
                end
            end
        end
    end

    % Charmatrix with names of filters
    mFilterNames = char(cFilter{:,1});
    mFilterFiles = char(cFilter{:,2});

    % Invoke the user interface
    hDialog = import_dialog(mFilterNames, mFilterFiles, sFilterDir);

    % Analyze Output
    if ~ishandle(hDialog)
        answer = 0;
    else
        handles = guidata(hDialog);
        answer = handles.answer;
        % OK pressed
        if answer == 1
            % Get the values from figure
            nFilter = get(handles.lstFilter, 'Value');
            % Remove figure from memory
            delete(hDialog);
            % Get the catalog file
            if isempty(FilePath)
                [newfile, newpath] = uigetfile('*.*', 'Choose catalog datafile');
            else
                [pathstr, name, ext] = fileparts(FilePath);
                newpath =[pathstr,filesep];
                newfile = [name,ext];
            end

            % Cancel pressed?
            if ~(isequal(newfile, 0) | isequal(newpath, 0))
                newfile = [newpath newfile];
                % Everything ok?
                if length(newfile) > 1
                    mCatalog = feval(cFilter{nFilter,2}, 1, newfile);
                end
            end
        else
            % Remove figure from memory
            delete(hDialog);
        end
    end

end

author: CelsoReyes


text:
function [mpc, warns] = psse_convert(warns, data, verbose)
%PSSE_CONVERT  Converts data read from PSS/E RAW file to MATPOWER case.
%   [MPC, WARNINGS] = PSSE_CONVERT(WARNINGS, DATA)
%   [MPC, WARNINGS] = PSSE_CONVERT(WARNINGS, DATA, VERBOSE)
%
%   Converts data read from a version RAW data file into a
%   MATPOWER case struct.
%
%   Input:
%       WARNINGS :  cell array of strings containing accumulated
%                   warning messages
%       DATA : struct read by PSSE_READ (see PSSE_READ for details).
%       VERBOSE :   1 to display progress info, 0 (default) otherwise
%
%   Output:
%       MPC : a MATPOWER case struct created from the PSS/E data
%       WARNINGS :  cell array of strings containing updated accumulated
%                   warning messages
%
%   See also PSSE_READ.

%   MATPOWER
%   Copyright (c) 2014-2016, Power Systems Engineering Research Center (PSERC)
%   by Yujia Zhu, PSERC ASU
%   and Ray Zimmerman, PSERC Cornell
%   Based on mpraw2mp.m, written by: Yujia Zhu, Jan 2014, yzhu54@asu.edu.
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% define named indices into bus, gen, branch 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;
[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;
[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;

%% options
sort_buses = 1;

%% defaults
if nargin < 3
    verbose = 0;
end
haveVlims = 0;
Vmin = 0.9;
Vmax = 1.1;

%%-----  case identification data  -----
baseMVA = data.id.SBASE;
rev = data.id.REV;

%%-----  bus data  -----
numbus = data.bus.num;
[nb, ncols] = size(numbus); %% number of buses, number of cols
bus = zeros(nb, VMIN);      %% initialize bus matrix
if rev < 24
    bus_name_col = 10;
else
    bus_name_col = 2;
end
if sort_buses
    [numbus, i] = sortrows(numbus, 1);
    bus_name = data.bus.txt(i, bus_name_col);
else
    bus_name = data.bus.txt(:, bus_name_col);
end
if rev < 24     %% includes loads
    bus(:, [BUS_I BUS_TYPE PD QD GS BS BUS_AREA VM VA BASE_KV ZONE]) = ...
        numbus(:, [1:9 11:12]);
elseif rev < 31     %% includes GL, BL
    bus(:, [BUS_I BASE_KV BUS_TYPE GS BS BUS_AREA ZONE VM VA]) = ...
        numbus(:, [1 3 4 5 6 7 8 9 10]);
else                %% fixed shunts and loads are in their own tables
    bus(:, [BUS_I BASE_KV BUS_TYPE BUS_AREA ZONE VM VA]) = ...
        numbus(:, [1 3 4 5 6 8 9]);
    if ncols >= 11 && all(all(~isnan(numbus(:, [10 11]))))
        haveVlims = 1;
        bus(:, [VMAX VMIN]) = numbus(:, [10 11]);
    end
end
if ~haveVlims  %% add default voltage magnitude limits if not provided
    warns{end+1} = sprintf('Using default voltage magnitude limits: VMIN = %g p.u., VMAX = %g p.u.', Vmin, Vmax);
    if verbose
        fprintf('WARNING: No bus voltage magnitude limits provided.\n         Using defaults: VMIN = %g p.u., VMAX = %g p.u.\n', Vmin, Vmax);
    end
    bus(:, VMIN) = Vmin;
    bus(:, VMAX) = Vmax;
end

%% create map of external bus numbers to bus indices
i2e = bus(:, BUS_I);
e2i = sparse(i2e, ones(nb, 1), 1:nb, max(i2e), 1);

%%-----  load data  -----
if rev >= 24
    nld = size(data.load.num, 1);
    loadbus = e2i(data.load.num(:,1));
    %% PSS/E loads are divided into:
    %%  1. constant MVA, (I=1)
    %%  2. constant current (I=I)
    %%  3. constant reactance/resistance (I = I^2)
    %% NOTE: reactive power component of constant admittance load is negative
    %%       quantity for inductive load and positive for capacitive load
    Pd = data.load.num(:,6) + data.load.num(:,8) .* bus(loadbus, VM) ...
            + data.load.num(:,10) .* bus(loadbus, VM).^2;
    Qd = data.load.num(:,7) + data.load.num(:,9) .* bus(loadbus, VM) ...
            - data.load.num(:,11) .* bus(loadbus, VM).^2;
    Cld = sparse(1:nld, loadbus, data.load.num(:,3), nld, nb);    %% only in-service-loads
    bus(:, [PD QD]) = Cld' * [Pd Qd];
end

%%-----  fixed shunt data  -----
if isfield(data, 'shunt')   %% rev > 30
    nsh = size(data.shunt.num, 1);
    shuntbus = e2i(data.shunt.num(:,1));
    Csh = sparse(1:nsh, shuntbus, data.shunt.num(:,3), nsh, nb);  %% only in-service shunts
    bus(:, [GS BS]) = Csh' * data.shunt.num(:, 4:5);
end

%%-----  switched shunt data  -----
nswsh = size(data.swshunt.num, 1);
swshuntbus = e2i(data.swshunt.num(:,1));
Cswsh = sparse(1:nswsh, swshuntbus, 1, nswsh, nb);
if rev <= 27
    bus(:, BS) = bus(:, BS) + Cswsh' * data.swshunt.num(:, 6);
elseif rev <= 29
    bus(:, BS) = bus(:, BS) + Cswsh' * data.swshunt.num(:, 7);
elseif rev < 32
    bus(:, BS) = bus(:, BS) + Cswsh' * data.swshunt.num(:, 8);
else
    bus(:, BS) = bus(:, BS) + Cswsh' * data.swshunt.num(:, 10);
end

%%-----  branch data  -----
nbr = size(data.branch.num, 1);
branch = zeros(nbr, ANGMAX);
branch(:, ANGMIN) = -360;
branch(:, ANGMAX) = 360;
branch(:, [F_BUS BR_R BR_X BR_B RATE_A RATE_B RATE_C]) = ...
    data.branch.num(:, [1 4 5 6 7 8 9]);
branch(:, T_BUS) = abs(data.branch.num(:, 2));  %% can be negative to indicate metered end
if rev <= 27        %% includes transformer ratio, angle
    branch(:, BR_STATUS) = data.branch.num(:, 16);
    branch(~isnan(data.branch.num(:, 10)), TAP) = ...
        data.branch.num(~isnan(data.branch.num(:, 10)), 10);
    branch(~isnan(data.branch.num(:, 11)), SHIFT) = ...
        data.branch.num(~isnan(data.branch.num(:, 11)), 11);
else
    branch(:, BR_STATUS)    = data.branch.num(:, 14);
end
%% integrate branch shunts (explicit shunts, not line-charging)
ibr = (1:nbr)';
fbus = e2i(branch(:, F_BUS));
tbus = e2i(branch(:, T_BUS));
nzf = find(fbus);               %% ignore branches with bad bus numbers
nzt = find(tbus);
if length(nzf) < nbr
    warns{end+1} = sprintf('%d branches have bad ''from'' bus numbers', nbr-length(nzf));
    if verbose
        fprintf('WARNING: %d branches have bad ''from'' bus numbers\n', nbr-length(nzf));
    end
end
if length(nzt) < nbr
    warns{end+1} = sprintf('%d branches have bad ''to'' bus numbers', nbr-length(nzt));
    if verbose
        fprintf('WARNING: %d branches have bad ''to'' bus numbers\n', nbr-length(nzt));
    end
end
Cf = sparse(ibr(nzf), fbus(nzf), branch(nzf, BR_STATUS), nbr, nb);  %% only in-service branches
Ct = sparse(ibr(nzt), tbus(nzt), branch(nzt, BR_STATUS), nbr, nb);  %% only in-service branches
if rev <= 27
    bus(:, [GS BS]) = bus(:, [GS BS]) + ...
        Cf' * data.branch.num(:, 12:13)*baseMVA + ...
        Ct' * data.branch.num(:, 14:15)*baseMVA;
else
    bus(:, [GS BS]) = bus(:, [GS BS]) + ...
        Cf' * data.branch.num(:, 10:11)*baseMVA + ...
        Ct' * data.branch.num(:, 12:13)*baseMVA;
end

%%-----  generator data  -----
ng = size(data.gen.num, 1);
genbus = e2i(data.gen.num(:,1));
gen = zeros(ng, APF);
gen(:, [GEN_BUS PG QG QMAX QMIN VG MBASE GEN_STATUS PMAX PMIN]) = ...
    data.gen.num(:, [1 3 4 5 6 7 9 15 17 18]);

%%-----  transformer data  -----
if rev > 27
    [transformer, bus, warns, bus_name] = psse_convert_xfmr(warns, data.trans2.num, data.trans3.num, verbose, baseMVA, bus, bus_name);
    branch = [branch; transformer];
end

%%-----  two-terminal DC transmission line data  -----
dcline = psse_convert_hvdc(data.twodc.num, bus);

%% assemble MPC
mpc = struct( ...
    'baseMVA',  baseMVA, ...
    'bus', bus, ...
    'bus_name', {bus_name}, ...
    'branch', branch, ...
    'gen', gen ...
);
if ~isempty(dcline)
    mpc.dcline = dcline;
    mpc = toggle_dcline(mpc, 'on');
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function xyl = Xl(s)
x = 0 ;
y = s ;

xyl = [x ; y] ;


text:
function P = spm_mesh_project(M, dat, method, varargin)
% Project volumetric data onto a mesh
% FORMAT P = spm_mesh_project(M, dat, method)
% M        - a patch structure, a handle to a patch
%            or a [nx3] vertices array
% dat      - a structure array [1xm] with fields dim, mat, XYZ and t
%            (see spm_render.m)
%            or a structure array [1xm] with fields mat and dat
%            or a structure array [1xm] from spm_vol.m
%            or a char array/cellstr of image filenames
% method   - interpolation method {'nn'}
% varargin - other parameters required by the interpolation method
%
% P        - a [mxn] projected data array
%__________________________________________________________________________
% Copyright (C) 2009-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_project.m 7543 2019-03-15 12:56:04Z guillaume $

if ishandle(M)
    V = get(M,'Vertices');
elseif isstruct(M) || isa(M,'gifti')
    V = M.vertices;
else
    V = M;
end

if nargin < 3, method = 'nn'; end
if ~strcmpi(method,'nn')
    error('Only Nearest Neighbours interpolation is available.');
end

if ischar(dat), dat = cellstr(dat); end
P = zeros(length(dat),size(V,1));
for i=1:numel(dat)
    if iscellstr(dat) || isfield(dat,'dt')
        if iscellstr(dat)
            v  = spm_vol(dat{i});
        else
            v  = dat(i);
        end
        Y      = spm_read_vols(v);
        mat    = v.mat;
    elseif isfield(dat,'dat')
        Y      = dat(i).dat;
        mat    = dat(i).mat;
    else
        if isfield(dat,'Z') %-xSPM structure
           dat = struct('dim',dat.DIM,'XYZ',dat.XYZ,'t',dat.Z,'mat',dat.M);
        end
        dim    = dat(i).dim(1:3);
        Y      = zeros(dim(:)');
        OFF    = dat(i).XYZ(1,:) + ...
                 dat(i).dim(1)*(dat(i).XYZ(2,:)-1 + ...
                 dat(i).dim(2)*(dat(i).XYZ(3,:)-1));
        Y(OFF) = dat(i).t .* (dat(i).t > 0);
        mat    = dat(i).mat;
    end
    XYZ        = double(inv(mat)*[V';ones(1,size(V,1))]);
    P(i,:)     = spm_sample_vol(Y,XYZ(1,:),XYZ(2,:),XYZ(3,:),0);
end


text:
function [anx,any,anz,dx,dy,dz,ierr] = focal_pt2nd(wpx,wpy,wpz,wtx,wty,wtz)

    % compute Cartesian component of P and T versors from outward normal and slip vectors
    %
    %     usage:
    %     call pt2nd(px,py,pz,tx,ty,tz,anx,any,anz,dx,dy,dz,ierr)
    %
    %     arguments:
    %     px,py,pz       components of P (maximum dilatation) axis vector
    %                    in the Aki-Richards Cartesian coordinate system (INPUT)
    %     tx,ty,tz       components of T (maximum tension) axis vector
    %                    in the Aki-Richards Cartesian coordinate system (INPUT)
    %     anx,any,anz    components of fault plane outward normal versor in the
    %                    Aki-Richards Cartesian coordinate system (OUTPUT)
    %     dx,dy,dz       components of slip versor in the Aki-Richards
    %                    Cartesian coordinate system (OUTPUT)
    %     ierr           error indicator (OUTPUT)
    %
    %     errors:
    %     1              input vectors not perpendicular among each other
    %
    %      call fpsset
    amistr=-360.;
    amastr=360.;
    amidip=0.;
    amadip=90.;
    amirak=-360.;
    amarak=360.;
    amitre=-360.;
    amatre=360.;
    amiplu=0.;
    amaplu=90.;
    orttol=2.;
    ovrtol=0.001;
    tentol=0.0001;
    dtor=0.017453292519943296;
    c360=360.;
    c90=90.;
    c0=0.;
    c1=1.;
    c2=2.;
    c3=3.;

    %%c
    anx=c0;
    any=c0;
    anz=c0;
    dx=c0;
    dy=c0;
    dz=c0;
    ierr=0;
    [ang] = focal_angle(wpx,wpy,wpz,wtx,wty,wtz);
    if (abs(ang - c90) > orttol)
        disp(['PT2ND: input vectors not perpendicular, angle=' num2str(ang)]);
        ierr = 1;
    end
    [pnorm,px,py,pz] = focal_norm(wpx,wpy,wpz);

    if (pz < c0)
        [px,py,pz] = focal_invert(px,py,pz);
    end
    [tnorm,tx,ty,tz] = focal_norm(wtx,wty,wtz);
    if (tz < c0)
        [tx,ty,tz] = focal_invert(tx,ty,tz);
    end
    anx = tx + px;
    any = ty + py;
    anz = tz + pz;
    [amn,anx,any,anz] = focal_norm(anx,any,anz);
    % c
    dx = tx - px;
    dy = ty - py;
    dz = tz - pz;
    [amn,dx,dy,dz] = focal_norm(dx,dy,dz);
    if(anz > c0)
        [anx,any,anz] = focal_invert(anx,any,anz);
        [dx,dy,dz] = focal_invert(dx,dy,dz);
    end
end


text:
function s = addlogi(s)
%ADDLOGI Add the logistic adistributions.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.9 $  $Date: 2004/01/24 09:35:06 $

j = length(s) + 1;
s(j).name = 'Logistic';
s(j).code = 'logistic';
s(j).pnames = {'mu' 'sigma'};
s(j).pdescription = {'location' 'scale'};
s(j).prequired = [false false];
s(j).fitfunc = @logifit;
s(j).likefunc = @logilike;
s(j).cdffunc = @logicdf;
s(j).pdffunc = @logipdf;
s(j).invfunc = @logiinv;
s(j).statfunc = @logistat;
s(j).loginvfunc = [];
s(j).logcdffunc = [];
s(j).hasconfbounds = false;
s(j).censoring = true;
s(j).paramvec = true;
s(j).support = [-Inf Inf];
s(j).closedbound = [false false];
s(j).iscontinuous = true;
s(j).islocscale = true;
s(j).uselogpp = false;

j = j + 1;
s(j).name = 'Log-Logistic';
s(j).code = 'loglogistic';
s(j).pnames = {'mu' 'sigma'};
s(j).pdescription = {'log location' 'log scale'};
s(j).prequired = [false false];
s(j).fitfunc = @loglfit;
s(j).likefunc = @logllike;
s(j).cdffunc = @loglcdf;
s(j).pdffunc = @loglpdf;
s(j).invfunc = @loglinv;
s(j).statfunc = @loglstat;
s(j).loginvfunc = @logiinv;
s(j).logcdffunc = @logicdf;
s(j).hasconfbounds = false;
s(j).censoring = true;
s(j).paramvec = true;
s(j).support = [0 Inf];
s(j).closedbound = [false false];
s(j).iscontinuous = true;
s(j).islocscale = true;
s(j).uselogpp = true;

% ==== Logistic distribution functions ====

% these distribution functions do not yet handle arrays of parameters

function y = logipdf(x, mu, sigma)
%LOGIPDF Logistic probability density function (pdf).
if (nargin<2), mu=0; end
if (nargin<3), sigma=1; end
sigma(sigma <= 0) = NaN;

z = (x - mu) ./ sigma;
k = (z>350); if any(k), z(k) = -z(k); end % prevent Inf/Inf
y = exp(z) ./ ((1 + exp(z)).^2 .* sigma);

function p = logicdf(x, mu, sigma)
%LOGICDF Logistic cumulative distribution function (cdf).
if (nargin<2), mu=0; end
if (nargin<3), sigma=1; end
sigma(sigma <= 0) = NaN;

p = 1 ./ (1 + exp(-(x - mu) ./ sigma));

function x = logiinv(p, mu, sigma)
%LOGIINV Inverse of the logistic cumulative distribution function (cdf).
if (nargin<2), mu=0; end
if (nargin<3), sigma=1; end
sigma(sigma <= 0) = NaN;

x = logit(p).*sigma + mu;

function r = logirnd(mu, sigma, varargin)
%LOGIRND Random arrays from the logistic distribution.
if (nargin<1), mu=0; end
if (nargin<2), sigma=1; end
sigma(sigma <= 0) = NaN;

[err, sizeOut] = statsizechk(2,mu,sigma,varargin{:});
if err > 0
    error('stats:logirnd:InconsistentSizes','Size information is inconsistent.');
end

p = rand(sizeOut);
r = log(p./(1-p)).*sigma + mu;

function [m,v] = logistat(mu, sigma)
%LOGISTAT Mean and variance for the logistic distribution.
if (nargin<1), mu=0; end
if (nargin<2), sigma=1; end
sigma(sigma <= 0) = NaN;

m = mu;
v = sigma.^2 .* pi.^2 ./ 3;

function [nlogL,acov] = logilike(params,data,cens,freq)
%LOGILIKE Negative log-likelihood for the logistic distribution.
if nargin < 4 || isempty(freq), freq = ones(size(data)); end
if nargin < 3 || isempty(cens), cens = zeros(size(data)); end

nlogL = logi_nloglf(params, data, cens, freq);
if nargout > 1
    acov = mlecov(params, data, 'nloglf',@logi_nloglf, 'cens',cens, 'freq',freq);
end

% ==== Logistic fitting functions ====

function [phat,pci] = logifit(x,alpha,cens,freq,opts)
%LOGIFIT Parameter estimates and confidence intervals for logistic data.

if nargin < 2 || isempty(alpha), alpha = .05; end
if nargin < 3 || isempty(cens), cens = zeros(size(x)); end
if nargin < 4 || isempty(freq), freq = ones(size(x)); end
if nargin < 5, opts = []; end

% Moment estimators as starting point
xunc = x(cens == 0);
start = [mean(xunc) std(xunc).*sqrt(3)./pi];

% The default options include turning statsfminbx's display off.  This
% function gives its own warning/error messages, and the caller can turn
% display on to get the text output from statsfminbx if desired.
options = statset(statset('logifit'), opts);
tolBnd = options.TolBnd;
options = optimset(options);
dfltOptions = struct('DerivativeCheck','off', 'HessMult',[], ...
    'HessPattern',ones(2,2), 'PrecondBandWidth',Inf, ...
    'TypicalX',ones(2,1), 'MaxPCGIter',1, 'TolPCG',0.1);

% Maximize the log-likelihood with respect to mu and sigma.
funfcn = {'fungrad' 'logifit' @logi_nloglf [] []};
[phat, nll, lagrange, err, output] = ...
         statsfminbx(funfcn, start, [-Inf; tolBnd], [Inf; Inf], ...
                     options, dfltOptions, 1, x, cens, freq);
if (err == 0)
    % statsfminbx may print its own output text; in any case give something
    % more statistical here, controllable via warning IDs.
    if output.funcCount >= options.MaxFunEvals
        wmsg = 'Maximum likelihood estimation did not converge.  Function evaluation limit exceeded.';
    else
        wmsg = 'Maximum likelihood estimation did not converge.  Iteration limit exceeded.';
    end
    warning('stats:logifit:IterOrEvalLimit',wmsg);
elseif (err < 0)
    error('stats:logifit:NoSolution',...
          'Unable to reach a maximum likelihood solution.');
end

if nargout > 1
    acov = mlecov(phat, x, 'nloglf',@logi_nloglf, 'cens',cens, 'freq',freq);
    probs = [alpha/2; 1-alpha/2];
    se = sqrt(diag(acov))';

    % Compute the CI for mu using a normal approximation for muhat.
    pci(:,1) = norminv(probs, phat(1), se(1));

    % Compute the CI for sigma using a normal approximation for
    % log(sigmahat), and transform back to the original scale.
    % se(log(sigmahat)) is se(sigmahat) / sigmahat.
    logsigci = norminv(probs, log(phat(2)), se(2)./phat(2));
    pci(:,2) = exp(logsigci);
end

function [nll,ngrad] = logi_nloglf(parms, x, cens, freq)
%LOGI_NLOGLF Objective function for logistic maximum likelihood.
mu = parms(1);
sigma = parms(2);
z = (x - mu) ./ sigma;
logitz = 1 ./ (1 + exp(-z));
clogitz = 1 ./ (1 + exp(z));
logclogitz = log(clogitz);
k = (z > 700); if any(k), logclogitz(k) = z(k); end % fix intermediate overflow

L = z + 2.*logclogitz - log(sigma);
ncen = sum(freq.*cens);
if ncen > 0
    cen = (cens == 1);
    L(cen) = logclogitz(cen);
end
nll = -sum(freq .* L);

if nargout > 1
    t = (2.*logitz - 1) ./ sigma;
    dL1 = t;
    dL2 = z.*t - 1./sigma;
    if ncen > 0
        t = logitz(cen) ./ sigma;
        dL1(cen) = t;
        dL2(cen) = z(cen) .* t;
    end
    ngrad = -[sum(freq .* dL1) sum(freq .* dL2)];
end

% ==== Log-Logistic distribution functions ====

% these distribution functions do not yet handle arrays of parameters

function y = loglpdf(x, mu, sigma)
%LOGLPDF Log-logistic probability density function (pdf).
if (nargin<2), mu=0; end
if (nargin<3), sigma=1; end
sigma(sigma <= 0) = NaN;

nonpos = (x <= 0);
x(nonpos) = realmin;
z = (log(x) - mu) ./ sigma;
c = ones(size(z));
k = (z>350); % prevent Inf/Inf
if any(k)
    z(k) = -z(k);
    c(k) = -1;
end
y = exp(z.*(1-c.*sigma) - mu) ./ ((1 + exp(z)).^2 .* sigma);
y(nonpos) = 0;
% the first and third of these would happen automatically for x==0, but
% generate LogOfZero warnings.  the second would be NaN.
y(x==0 & sigma<1) = 0;
y(x==0 & sigma==1) = 1;
y(x==0 & sigma>1) = Inf;

function p = loglcdf(x, mu, sigma)
%LOGLCDF Log-logistic cumulative distribution function (cdf).
if (nargin<2), mu=0; end
if (nargin<3), sigma=1; end
sigma(sigma <= 0) = NaN;

nonpos = (x <= 0);
x(nonpos) = realmin;
p = 1 ./ (1 + exp(-(log(x) - mu) ./ sigma));
% this would happen automatically for x==0, but generates LogOfZero warnings
p(nonpos) = 0;

function x = loglinv(p, mu, sigma)
%LOGLINV Inverse of the log-logistic cumulative distribution function (cdf).
if (nargin<2), mu=0; end
if (nargin<3), sigma=1; end
sigma(sigma <= 0) = NaN;

x = exp(logit(p).*sigma + mu);

function r = loglrnd(mu, sigma, varargin)
%LOGLRND Random arrays from the log-logistic distribution.
if (nargin<1), mu=0; end
if (nargin<2), sigma=1; end
sigma(sigma <= 0) = NaN;

[err, sizeOut] = statsizechk(2,mu,sigma,varargin{:});
if err > 0
    error('stats:loglrnd:InconsistentSizes','Size information is inconsistent.');
end

p = rand(sizeOut);
r = exp(log(p./(1-p)).*sigma + mu);

function [m,v] = loglstat(mu, sigma)
%LOGLSTAT Mean and variance for the log-logistic distribution.
if (nargin<1), mu=0; end
if (nargin<2), sigma=1; end
sigma(sigma <= 0) = NaN;

if sigma < 1
    m = exp(mu + gammaln(1+sigma) + gammaln(1-sigma));
else
    m = Inf;
end
if sigma < .5
    v = exp(2.*mu + gammaln(1+2.*sigma) + gammaln(1-2.*sigma)) - m.^2;
else
    v = Inf;
end

function [nlogL,acov] = logllike(params,data,cens,freq)
%LOGLLIKE Negative log-likelihood for the log-logistic distribution.
if nargin < 4 || isempty(freq), freq = ones(size(data)); end
if nargin < 3 || isempty(cens), cens = zeros(size(data)); end

nlogL = logl_nloglf(params, data, cens, freq);
if nargout > 1
    acov = mlecov(params, data, 'nloglf',@logl_nloglf, 'cens',cens, 'freq',freq);
end

% ==== Log-Logistic fitting functions ====

function [phat,pci] = loglfit(x,alpha,cens,freq,opts)
%LOGLFIT Parameter estimates and confidence intervals for log-logistic data.

if nargin < 2 || isempty(alpha), alpha = .05; end
if nargin < 3 || isempty(cens), cens = zeros(size(x)); end
if nargin < 4 || isempty(freq), freq = ones(size(x)); end
if nargin < 5, opts = []; end

if any(x <= 0)
    error('stats:loglfit:BadData','The data in X must be positive');
end

% Moment estimators as starting point
logxunc = log(x(cens == 0));
start = [mean(logxunc) std(logxunc).*sqrt(3)./pi];

% The default options include turning statsfminbx's display off.  This
% function gives its own warning/error messages, and the caller can turn
% display on to get the text output from statsfminbx if desired.
options = statset(statset('loglfit'), opts);
tolBnd = options.TolBnd;
options = optimset(options);
dfltOptions = struct('DerivativeCheck','off', 'HessMult',[], ...
    'HessPattern',ones(2,2), 'PrecondBandWidth',Inf, ...
    'TypicalX',ones(2,1), 'MaxPCGIter',1, 'TolPCG',0.1);

% Maximize the log-likelihood with respect to mu and sigma.
funfcn = {'fungrad' 'loglfit' @logl_nloglf [] []};
[phat, nll, lagrange, err, output] = ...
         statsfminbx(funfcn, start, [-Inf; tolBnd], [Inf; Inf], ...
                     options, dfltOptions, 1, x, cens, freq);
if (err == 0)
    % statsfminbx may print its own output text; in any case give something
    % more statistical here, controllable via warning IDs.
    if output.funcCount >= options.MaxFunEvals
        wmsg = 'Maximum likelihood estimation did not converge.  Function evaluation limit exceeded.';
    else
        wmsg = 'Maximum likelihood estimation did not converge.  Iteration limit exceeded.';
    end
    warning('stats:loglfit:IterOrEvalLimit',wmsg);
elseif (err < 0)
    error('stats:loglfit:NoSolution',...
          'Unable to reach a maximum likelihood solution.');
end

if nargout > 1
    acov = mlecov(phat, x, 'nloglf',@logi_nloglf, 'cens',cens, 'freq',freq);
    probs = [alpha/2; 1-alpha/2];
    se = sqrt(diag(acov))';

    % Compute the CI for mu using a normal approximation for muhat.
    pci(:,1) = norminv(probs, phat(1), se(1));

    % Compute the CI for sigma using a normal approximation for
    % log(sigmahat), and transform back to the original scale.
    % se(log(sigmahat)) is se(sigmahat) / sigmahat.
    logsigci = norminv(probs, log(phat(2)), se(2)./phat(2));
    pci(:,2) = exp(logsigci);
end

function [nll,ngrad] = logl_nloglf(parms, x, cens, freq)
%LOGL_NLOGLF Objective function for log-logistic maximum likelihood.
mu = parms(1);
sigma = parms(2);
logx = log(x);
z = (logx - mu) ./ sigma;
logitz = 1 ./ (1 + exp(-z));
clogitz = 1 ./ (1 + exp(z));
logclogitz = log(clogitz);
k = (z > 700); if any(k), logclogitz(k) = z(k); end % fix intermediate overflow

L = z + 2.*logclogitz - log(sigma) - logx;
ncen = sum(freq.*cens);
if ncen > 0
    cen = (cens == 1);
    L(cen) = logclogitz(cen);
end
nll = -sum(freq .* L);

if nargout > 1
    t = (2.*logitz - 1) ./ sigma;
    dL1 = t;
    dL2 = z.*t - 1./sigma;
    if ncen > 0
        t = logitz(cen) ./ sigma;
        dL1(cen) = t;
        dL2(cen) = z(cen) .* t;
    end
    ngrad = -[sum(freq .* dL1) sum(freq .* dL2)];
end

% ==== utility functions ====

function logitp = logit(p)
%LOGIT Logistic transformation, handling edge and out of range.
logitp = repmat(NaN,size(p));
logitp(p==0) = -Inf;
logitp(p==1) = Inf;
ok = (0<p & p<1);
logitp(ok) = log(p(ok)./(1-p(ok)));

author: zouchuhang

repo: LayoutNet