text:
function varargout = plot(f, varargin)
%PLOT   Plot a BALLFUN on the ball
%   PLOT(F) creates a plot showing F on the ball.
%
%   PLOT(f, 'WedgeAz') plot a BALLFUN with a wedge in the azimuthal
%   (longitude) direction removed.
%
%   PLOT(f, 'WedgePol') plot a BALLFUN with a wedge in the polar
%   (latitude) direction removed.
%
% EXAMPLES:
%   f = cheb.galleryball;
%   plot(f)
%   plot(f, 'WedgeAz')
%   plot(f, 'WedgePol')
%
% See also BALLFUN/SURF, BALLFUN/SLICE.

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

% Check if the function is empty
if isempty(f)
    error('CHEBFUN:BALLFUN:plot:isempty','Function is empty.');
end

% Add a warning of the function is not real
if (f.isReal == 0)
    warning('CHEBFUN:BALLFUN:plot:isReal','Function is not real, plotting the real part.');
end

if ( nargin == 1 )
    h = plotBall(f);
elseif ( nargin == 2 ) && ( strcmpi(varargin{1},'wedgeaz') )
    h = plotWedgeAz(f);
elseif ( nargin == 2 ) && ( strcmpi(varargin{1},'wedgepol') )
    h = plotWedgePol(f);
else
    error('CHEBFUN:BALLFUN:plot:input Invalid input arguments')
end

if ( nargout > 0 )
    varargout = { h };
end
end

function h = plotBall(f)
% Plot a BALLFUN function on the ball

% Is the plot currently being held?
plotOnHold = ishold;
% Default plotting options
defaultOpts = {'facecolor', 'interp','edgecolor', 'none'};

% Define the size of F:
[m,n,p] = size(f);

% m >= 25 and n, p >= 28
m = 25*(m < 25) + m*(m >= 25);
n = 28*(n < 28) + n*(n>=28);
p = 28*(p < 28) + p*(p>=28);

% Impose m = 1 [6] and n, p = 0 [4] to avoid errors in the plot
m = m + mod(1-mod(m,6),6);
n = n + mod(4-mod(n,4),4);
p = p + mod(4-mod(p,4),4);

% Get the coeffs of the ballfun function F(r,lam,th), lam is the
% azimuthal variable in [-pi,pi] and theta the polar variable in [0,pi]
F = coeffs3(f, m, n, p);

% Convert to values
ff = real(ballfun.coeffs2vals(F));

% Permute lambda and theta
ff = permute(ff,[1 3 2]);

% Evaluation points
r   = chebpts( m );
lam  = [pi*trigpts( n ); pi];
th = [pi*trigpts( p );pi]-pi/2;

% Remove doubled-up data
r = r(floor(m/2)+1:end);
th = th(floor(p/2)+1:end);

% Reverse theta : 1st element of the array is theta = pi (South Pole), last element is
% th = 0 (not included) (North Pole)
ff  = ff(floor(m/2)+1:end,[1 end:-1:floor(p/2)+1],:);
ff(:,:,end+1) = ff(:,:,1);

% Define the meshgrid
[tt, rr, ll] = meshgrid(th, r, lam);

% Slices in the cylinder to plot
% Find the indice of r = 0.5
[~,idr]=min(abs(r-0.5));
rslice = rr(idr,1,1);
tslice = tt(1,[1,floor(p/4)+1],1);
lslice = ll(1,1,[1,floor(n/4)+1]);

hslicer = slice(tt,rr,ll,ff,tslice,rslice,lslice);

hold on
for j = 1:numel(hslicer)
    h = hslicer(j);
    [xs,ys,zs] = sph2cart(h.ZData,h.XData,h.YData);
    h = surf(xs,ys,zs,h.CData,defaultOpts{:});
end
delete(hslicer);

if ~plotOnHold
    hold off;
end

axis([-1 1 -1 1 -1 1])
daspect([1 1 1])

camlight('headlight');
lighting phong;
material dull;
end

function h = plotWedgeAz(f)
% Plot a BALLFUN with a wedge in the azimuthal direction removed

% Is the plot currently being held?
plotOnHold = ishold;
% Default plotting options
defaultOpts = {'facecolor', 'interp','edgecolor', 'none'};

% Azimuthal (longitude) values to include.  TODO: Make these optional inputs
az_intvl = [-pi/2 pi];

% Define the size of f:
[m,n,p] = size(f);

% m >= 25 and n, p >= 28
m = 25*(m < 25) + m*(m >= 25);
n = 28*(n < 28) + n*(n>=28);
p = 28*(p < 28) + p*(p>=28);

% Impose m = 1 [6] and n, p = 0 [4] to avoid errors in the plot
m = m + mod(1-mod(m,6),6);
n = n + mod(4-mod(n,4),4);
p = p + mod(4-mod(p,4),4);

% Construct the values of lambda and theta to plot on the outer sphere (r=1)
lam = linspace(az_intvl(1),az_intvl(2),n);
th = linspace(0,pi,p)';

% Evaluate the function on the outer sphere
ff = permute(fevalm(f,1,lam,th),[3 2 1]);

% Plot the result
h = surf(sin(th)*cos(lam),sin(th)*sin(lam),cos(th)*ones(1,n),ff,defaultOpts{:});
hold on

% Construct the values of r and theta to plot from the origin to the outer
% sphere (r=1).
r = chebpts(m); r = r(floor(m/2)+1:end); th = th';

% Evaluate the function on the wedge from the origin along the az_intvl(1).
ff = permute(fevalm(f,r,lam(1),th),[1 3 2]);
% Plot the result
surf(r*sin(th)*cos(lam(1)),r*sin(th)*sin(lam(1)),r*cos(th),ff,defaultOpts{:})

% Evaluate the function on the wedge from the origin along the az_intvl(2).
ff = permute(fevalm(f,r,lam(end),th),[1 3 2]);
% Plot the result
surf(r*sin(th)*cos(lam(end)),r*sin(th)*sin(lam(end)),r*cos(th),ff,defaultOpts{:})

if ~plotOnHold
    hold off;
end

camlight('headlight');
lighting phong;
material dull;

axis([-1 1 -1 1 -1 1])
daspect([1 1 1])
end

function h = plotWedgePol(f)
% Plot a BALLFUN with a wedge in the polar direction removed

% Is the plot currently being held?
plotOnHold = ishold;
% Default plotting options
defaultOpts = {'facecolor', 'interp','edgecolor', 'none'};

% Polar (latitude) values to include.  TODO: Make these optional inputs
pol_intvl = [pi/2 pi];
% Azimuthal (longitude) values to include.  TODO: Make these optional inputs
az_intvl = [0 pi];

% Define the size of f:
[m,n,p] = size(f);

% m >= 25 and n, p >= 28
m = 25*(m < 25) + m*(m >= 25);
n = 28*(n < 28) + n*(n>=28);
p = 28*(p < 28) + p*(p>=28);

% Impose m = 1 [6] and n, p = 0 [4] to avoid errors in the plot
m = m + mod(1-mod(m,6),6);
n = n + mod(4-mod(n,4),4);
p = p + mod(4-mod(p,4),4);

% Construct the values of lambda and theta to plot on the outer sphere (r=1)
% where the sphere is closed
lam = linspace(-pi,pi,n);
th = linspace(pol_intvl(1),pol_intvl(2),p)';

% Evaluate the function on the outer sphere
ff = permute(fevalm(f,1,lam,th),[3 2 1]);

% Plot the result
h = surf(sin(th)*cos(lam),sin(th)*sin(lam),cos(th)*ones(1,n),ff,defaultOpts{:});
hold on

% Construct the values of lambda and theta to plot on the outer sphere (r=1)
% where the sphere is open at the wedge
lam = linspace(az_intvl(1),az_intvl(2),n);
th = linspace(0,pol_intvl(1),p)';

% Evaluate the function on the outer sphere
ff = permute(fevalm(f,1,lam,th),[3 2 1]);

% Plot the result
surf(sin(th)*cos(lam),sin(th)*sin(lam),cos(th)*ones(1,n),ff,defaultOpts{:})

% Construct the values of r and theta to plot from the origin to the outer
% sphere (r=1).
r = chebpts(m); r = r(floor(m/2)+1:end); th = linspace(0,pol_intvl(1),p);

% Evaluate the function on the wedge from the origin along the az_intvl(1).
ff = permute(fevalm(f,r,lam(1),th),[1 3 2]);
% Plot the result
surf(r*sin(th)*cos(lam(1)),r*sin(th)*sin(lam(1)),r*cos(th),ff,defaultOpts{:})

% Construct the values of r and theta to plot from the origin to the outer
% sphere (r=1).
r = chebpts(m); r = r(floor(m/2)+1:end); th = linspace(0,pol_intvl(1),p);

% Evaluate the function on the wedge from the origin along the az_intvl(1).
ff = permute(fevalm(f,r,lam(end),th),[1 3 2]);
% Plot the result
surf(r*sin(th)*cos(lam(end)),r*sin(th)*sin(lam(end)),r*cos(th),ff,defaultOpts{:})

% Construct the values of r and lambda to plot from the origin to the outer
% sphere (r=1): slice along
r = chebpts(m); r = r(floor(m/2)+1:end); lam = linspace(-pi,pi,p);

% Evaluate the function on the wedge from the origin along the az_intvl(2).
ff = fevalm(f,r,lam,pol_intvl(1));
% Plot the result
surf(r*sin(pol_intvl(1))*cos(lam),r*sin(pol_intvl(1))*sin(lam),r*cos(pol_intvl(1))*ones(size(lam)),ff,defaultOpts{:})

if ~plotOnHold
    hold off;
end

camlight('headlight');
lighting phong;
material dull;

axis([-1 1 -1 1 -1 1])
daspect([1 1 1])
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% example_1 - inversion of a irrational fraction F(s)
clear, close all
[t1,ft1]=INVLAP('tanh(s)/s',0.01,20,1000);
[t2,ft2]=INVLAP('tanh(s)/s',0.01,20,2000,6,280,59);
figure(3)
set(3,'color','white')
subplot(2,1,1)
plot(t1,ft1), grid on, zoom on
xlabel('t [s]'), ylabel('f(t)')
title('rectangular periodic wave')
subplot(2,1,2)
plot(t2,ft2), grid on, zoom on
xlabel('t [s]'), ylabel('f(t)')
title('improved accuracy')

author: Sable

repo: mcbench-benchmarks


text:
close all
clear all
clc

x = -pi:.1:pi;
y = sin(x);
p = plot(x,y)
set(gca,'XTick',-pi:pi/2:pi)
set(gca,'XTickLabel',{'-pi','-pi/2','0','pi/2','pi'})
xlabel('-\pi \leq \Theta \leq \pi')
ylabel('sin(\Theta)')
title('Simulation Results')
text(-pi/4,sin(-pi/4),'\leftarrow sin(-\pi\div4)',...
     'HorizontalAlignment','left')
set(p,'Color','red','LineWidth',2)

latex_fig(10, 3, 1.5)


text:
function node_xyz = sphere_grid_q16_node_xyz ( nelemx, nelemy )

%*****************************************************************************80
%
%% SPHERE_GRID_Q16_NODE_XYZ produces node coordinates for a Q16 sphere grid.
%
%  Discussion:
%
%    The number of nodes to be generated is
%
%      NODE_NUM = 9 * NELEMX * NELEMY - 3 * NELEMX + 2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.
%
%    Output, real NODE_XYZ(3,NODE_NUM), the node coordinates.
%
  node = 0;

  for j = 3 * nelemy + 1 : -1 : 1

    phi = ( j - 1 ) * pi / ( 3 * nelemy );

    if ( j == 1 )

      node = node + 1;
      node_xyz(1,node) =  0.0;
      node_xyz(2,node) =  0.0;
      node_xyz(3,node) =  1.0;

    elseif ( j < 3 * nelemy + 1 )

      for i = 1 : 3 * nelemx

        theta = ( i - 1 ) * 2.0 * pi / ( 3 * nelemx );

        node = node + 1;
        node_xyz(1,node) = cos ( theta ) * sin ( phi );
        node_xyz(2,node) = sin ( theta ) * sin ( phi );
        node_xyz(3,node) =                 cos ( phi );

      end

    elseif ( j == 3 * nelemy + 1 )

      node = node + 1;
      node_xyz(1,node) =  0.0;
      node_xyz(2,node) =  0.0;
      node_xyz(3,node) = -1.0;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ n_data, x, fx ] = fresnel_sin_values ( n_data )

%*****************************************************************************80
%
%% FRESNEL_SIN_VALUES returns some values of the Fresnel sine integral function.
%
%  Discussion:
%
%    The Fresnel sine integral is defined by
%
%      S(X) = integral ( 0 <= T <= X ) sin ( pi * T^2 / 2 ) / T dT
%
%    In Mathematica, the function can be evaluated by:
%
%      FresnelS[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 16;

  fx_vec = [ ...
     0.0000000000000000E+00, ...
     0.4187609161656762E-02, ...
     0.3335943266061318E-01, ...
     0.1105402073593870E+00, ...
     0.2493413930539178E+00, ...
     0.4382591473903548E+00, ...
     0.6234009185462497E+00, ...
     0.7135250773634121E+00, ...
     0.6388876835093809E+00, ...
     0.4509387692675831E+00, ...
     0.3434156783636982E+00, ...
     0.4557046121246569E+00, ...
     0.6196899649456836E+00, ...
     0.5499893231527195E+00, ...
     0.3915284435431718E+00, ...
     0.4963129989673750E+00 ];

  x_vec = [ ...
     0.0E+00, ...
     0.2E+00, ...
     0.4E+00, ...
     0.6E+00, ...
     0.8E+00, ...
     1.0E+00, ...
     1.2E+00, ...
     1.4E+00, ...
     1.6E+00, ...
     1.8E+00, ...
     2.0E+00, ...
     2.2E+00, ...
     2.4E+00, ...
     2.6E+00, ...
     2.8E+00, ...
     3.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

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

  return
end


text:
% Upsampling procedure.
%
% Argments:
%   'I': greyscale image
%   'odd': 2-vector of binary values, indicates whether the upsampled image
%   should have odd size for the respective dimensions
%   'filter': upsampling filter
%
% If image width W is odd, then the resulting image will have width (W-1)/2+1,
% Same for height.
%
% tom.mertens@gmail.com, August 2007
%

function R = upsample_(I,odd,filter)

% increase resolution
I = padarray(I,[1 1 0],'replicate'); % pad the image with a 1-pixel border
r = 2*size(I,1);
c = 2*size(I,2);
k = size(I,3);
R = zeros(r,c,k);
R(1:2:r, 1:2:c, :) = 4*I; % increase size 2 times; the padding is now 2 pixels wide

% interpolate, convolve with separable filter
R = imfilter(R,filter);     %horizontal
R = imfilter(R,filter');    %vertical

% remove the border
R = R(3:r - 2 - odd(1), 3:c - 2 - odd(2), :);

author: mahmoudnafifi


text:
% NCUGRID Extention of dataset object class for unstructured grid datasets.
% Usage:
%            >> dataset = ncugrid(uri);
%
% NCTOOLBOX (http://code.google.com/p/nctoolbox)
classdef ncugrid < handle

    properties (SetAccess = private)
        dataset % ncdataset obj (may not even be needed) since I am reproducing essential functions in
        % this ncugrid class...
        netcdfugrid % Kyle's netcdf u grid dataset
        numberofgrids % total number of meshes
        meshes % list of mesh objs
        cells % number of cells
        nodes % number of nodes
        edges % number of edges
        faces % number of faces
        variables % list of variable names
    end

    methods

        function obj = ncugrid(nc)
            echo off
            %             import java.io.IOException;
            import java.util.Formatter;
            %             import java.util.List;
            % %             import ucar.*
            %             import ucar.nc2.Variable;
            import ucar.nc2.constants.FeatureType;
            %             import ucar.nc2.dt.UGridDataset;
            %             import ucar.nc2.dt.ugrid.Cell;
            %             import ucar.nc2.dt.ugrid.Entity;
            %             import ucar.nc2.dt.ugrid.Mesh;
            % %             import ucar.nc2.dt.ugrid.topology.*;
            %             import ucar.nc2.dt.ugrid.geom.LatLonPoint2D;
            import ucar.nc2.ft.FeatureDatasetFactoryManager;
            import ucar.nc2.util.CancelTask;
            %             import ucar.unidata.geoloc.LatLonPoint;
            %             import ucar.unidata.geoloc.LatLonPointImpl;
            %             import ucar.unidata.geoloc.LatLonRect;
            import java.io.IOException;
            import java.util.ArrayList;
            import java.util.List;
            import ucar.ma2.Array;
            import ucar.ma2.DataType;
            import ucar.nc2.Attribute;
            import ucar.nc2.Dimension;
            import ucar.nc2.Variable;
            import ucar.nc2.dataset.CoordinateSystem;
            import ucar.nc2.dataset.NetcdfDataset;
            import ucar.nc2.dataset.VariableDS;
            import ucar.nc2.dt.ugrid.Cell;
            import ucar.nc2.dt.ugrid.Edge;
            import ucar.nc2.dt.ugrid.Face;
            import ucar.nc2.dt.ugrid.Node;
            import ucar.nc2.dt.ugrid.UGridDataset;
            import ucar.nc2.dt.ugrid.geom.LatLonPoint2D;

            if ischar(nc)
                try
                    obj.dataset = ncdataset(nc);  % src is a string URL/File
                catch % if it doesnt work maybe try dods instead of http
                    nc(1:4) = 'dods';
                    obj.dataset = ncdataset(nc);  % src is a string URL/File
                end
                cancelTask = [];
                %                 try
                obj.netcdfugrid = FeatureDatasetFactoryManager.open(FeatureType.UGRID, nc, cancelTask, Formatter());
                %                 catch me
                %                     me.throw
                %                 end

            elseif isa(nc, 'ncdataset')
                obj.dataset = nc;             % src is a ncdataset
                cancelTask = [];
                obj.netcdfugrid = FeatureDatasetFactoryManager.open(FeatureType.UGRID, nc.link, cancelTask, Formatter());
            elseif isa(nc,'ucar.nc2.dt.ugrid.UGridDataset')
                obj.dataset = nc;             % src is a raw ugriddataset in memory
                cancelTask = [];
                obj.netcdfugrid = FeatureDatasetFactoryManager.open(FeatureType.UGRID, nc, cancelTask, Formatter());
            else
                ex = MException('NCUGRID:ncugrid', 'Invalid dataset was specified');
                ex.throw;
            end

            %set properties
            u = obj.netcdfugrid;

%             vars = u.getMeshVariables();
            vars = obj.dataset.variables;
%             n = size(vars);
%             obj.variables = cell(n, 1);
%             for i = 1:(n)
%                 obj.variables{i, 1} = char(vars.get(i - 1).getName());
%             end
           obj.variables = vars;

            meshsets = u.getMeshsets();
            %       UGridDataset.getMeshsets().get(0).getMesh()
            obj.meshes = meshsets;
            obj.numberofgrids = meshsets.size;
            for i = 1:(obj.numberofgrids)
                obj.cells(i, 1) = meshsets.get(i - 1).getMesh.getSize();
                obj.nodes(i, 1) = meshsets.get(i - 1).getMesh.getNodeSize();
                obj.edges(i, 1) = meshsets.get(i - 1).getMesh.getEdgeSize();
                obj.faces(i, 1) = meshsets.get(i - 1).getMesh.getFaceSize();
            end

        end % end constructor

        function uvar = uvariable(obj, varName)
            % NCUGRID.UVARIABLE - Instantiates a ugrid variable object from the variable 'varName'.
            % Usage:
            %             uvar_obj = nc.uvariable(varName);
            %
            uvar = ncuvariable(obj, varName);
        end % uvariable end

        function ss = variableSubset(obj, varName, struct) % input subset structure
                    % NCUGRID.unstructuredLatLonSubset  - Function to subset an unstructured model grid by lat/lon
                    % bounding box using subsetting methods in the ugrid-java.
                    %
                    % Usage:
                    %              d = nc.unstructuredLatLonSubset(variableName, subsetstructure)
                    %
                    % Subset-structure:
                    %              subsetstructure.lat = [minlat maxlat];
                    %              subsetstructure.lon = [minlon maxlon];
                    %
                    import ucar.nc2.dt.ugrid.geom.LatLonPolygon2D;
                    import ucar.nc2.dt.ugrid.geom.LatLonRectangle2D;
                    import ucar.unidata.geoloc.LatLonPoint;
                    import ucar.unidata.geoloc.LatLonPointImpl;
                    import ucar.unidata.geoloc.LatLonRect;
%                     if nargin ==3
%                         struct = varargin{2};
%                         varName = varargin{1};
%                         % Find the mesh that contains the variable of interest.
%                         mesh_ident = obj.attribute('mesh', varName);
%                         mesh_vars = regexp(mesh_ident, ' ', 'split');
%                         mesh_ind = str2double(mesh_vars{1}(5:end));
%                     elseif nargin == 2
%                         struct = varargin{1};
%                     end
                    % Need to construct lat/lon rect:
                    maxlat  = max(struct.lat);
                    minlon = min(struct.lon);
                    minlat = min(struct.lat);
                    maxlon = max(struct.lon);
                    bbox = LatLonRect(LatLonPointImpl(maxlat, minlon), LatLonPointImpl(minlat, maxlon));

%                     if nargin == 3
%                         % Get all the meshes, access the correct mesh for the variable, then call subset method on mesh.
%                         meshes = obj.meshes;
%                         mesh = meshes.get(mesh_ind-1).getMesh();
%                         if mesh.getTreeSize() == 0
%                             mesh.buildRTree();
%                         end
%         %                 v = obj.netcdfugrid.getMeshVariableByName(varName);
%         %                 subsat = v.subsetToDataset(bbox); % Subsat is ugriddataset
%                         v = obj.netcdfugrid;
%                         subsat = v.subset(bbox); % Subsat is ugriddataset
%                     elseif nargin == 2
%                         v = obj.netcdfugrid;
%                         subsat = v.subset(bbox); % Subsat is ugriddataset
%                     end
                    mesh  = obj.netcdfugrid.getMeshVariableByName(varName).getMeshset().getMesh();
%                     mesh = meshvariable.getMeshset.getMesh();
                    if mesh.getTreeSize() == 0
                        mesh.buildRTree();
                    else
                        % pass
                    end
%                     clear meshvariable
%                     meshvariable = obj.netcdfugrid.getMeshVariableByName(varName);
                    subsatdataset = obj.netcdfugrid.getMeshVariableByName(varName).subsetToSelf(bbox);
%                     disp(subsatdataset.getNetcdfDataset.toString())
                    % Get netcdfj variable class and read all (since we already subset)
                    submeshv = subsatdataset.getMeshVariableByName(varName);
                    connect_meshset = submeshv.getMeshset();
%                     temp = connect_meshset.getMesh();
%                     temp.toString()

                    % Get variable values from subsat meshVariable
                    submeshv = submeshv.getVariable();
                    array = submeshv.read();
                    values = array.copyToNDJavaArray();
                    ss.data = values;
                    % Related coordinates and connectivity
                    ss.grid = obj.getSomeCoordinateData(varName, subsatdataset, connect_meshset);
        %         end

        end % end subset

        function sub = geosubset(obj, struct)
            import ucar.nc2.dt.ugrid.geom.LatLonPolygon2D;
            import ucar.nc2.dt.ugrid.geom.LatLonRectangle2D;
            import ucar.unidata.geoloc.LatLonPoint;
            import ucar.unidata.geoloc.LatLonPointImpl;
            import ucar.unidata.geoloc.LatLonRect;
            %                     if nargin ==3
            %                         struct = varargin{2};
            %                         varName = varargin{1};
            %                         % Find the mesh that contains the variable of interest.
            %                         mesh_ident = obj.attribute('mesh', varName);
            %                         mesh_vars = regexp(mesh_ident, ' ', 'split');
            %                         mesh_ind = str2double(mesh_vars{1}(5:end));
            %                     elseif nargin == 2
            %                         struct = varargin{1};
            %                     end
            % Need to construct lat/lon rect:
            maxlat  = max(struct.lat);
            minlon = min(struct.lon);
            minlat = min(struct.lat);
            maxlon = max(struct.lon);
            bbox = LatLonRect(LatLonPointImpl(maxlat, minlon), LatLonPointImpl(minlat, maxlon));
            subsatdataset = obj.netcdfugrid.subset(bbox);
%             disp(subsatdataset.getNetcdfDataset.toString())
        end % end geosubset (datasetSubset)

        function a = attributes(obj, variable)
            % NCUGRID.ATTRIBUTES returns the attributes of the variable as an
            % n x 2 cell array.
            %
            % Use as:
            %   a = ncdataset.attributes
            %   a = ncdataset.attributes(variableName)
            %
            % Inputs:
            %   variableName = The name of the variable whos attributes you want
            %       to retrieve. If no argument is specified then the
            %       global attributes are returned.
            %
            % Return:
            %   An (n, 2) cell array. The 1st column contains the attribute name
            %   The 2nd column contains the attribute value. Each row is essentially
            %   a key-value pair.
            %
            % Hints:
            %   Here's how to obtain cell arrays of the keys and corresponding values
            %   and search for a particular key (or attribute name)
            %     ds = ncugrid('http://somewhere/data.nc');
            %     attr = ds.attributes('time');
            %     keys = attr(:, 1);    % Cell array of keys
            %     values = attr(:, 2);  % Cell array of values
            %     i = find(ismember(keys, 'units')); % search for units attribute
            %     units = values{i};  % Retrieve the units value
            %             if (nargin < 2)
            %                 % Get global attributes
            %                 aa = obj.netcdfugrid.getGlobalAttributes();
            %             else
            %                 % Get attributes for the variable
            %                 variable = ucar.nc2.NetcdfFile.escapeName(variable);
            %                 v = obj.netcdfugrid.getNetcdfFile().findVariable(variable);
            %                 if isempty(v)
            %                     warning('NCTOOLBOX:ncugrid:attributes', ['Could not find the variable ']);
            %                 end
            %                 aa = v.getAttributes();
            %             end
            %
            %             if ~isempty(aa)
            %                 n = aa.size();
            %                 a = cell(n, 2);
            %                 for i = 1:n
            %                     at = aa.get(i - 1);
            %                     a{i, 1} = char(at.getName());
            %                     if (at.isString())
            %                         a{i, 2} = char(at.getStringValue());
            %                     else
            %                         a{i, 2} = at.getValues().copyToNDJavaArray();
            %                     end
            %                 end
            %             else
            %                 % Show warning, return empty cell array
            %                 warning('NCTOOLBOX:ncugrid:attributes', 'No attributes were found');
            %                 a = cell(1, 2);
            %             end
            if nargin == 2
                a = obj.dataset.attributes(variable);
            elseif nargin == 1

                a = obj.dataset.attributes();
            else
                % pass
            end
        end % end attributes

        function val = attribute(obj, varargin)
            % NCUGRID.ATTRIBUTE returns the value a global attribute specified by its key or the
            % variable attribute specified by key and variable.
            %
            % Use as:
            %   a = ncdataset.attribute('title')
            %   a = ncdataset.attribute(key)
            %
            %   a = ncdataset.attribute('temp', 'title')
            %   a = ncdataset.attribute(variableName, key)
            %
            % Inputs:
            %   key = The name of the attribute field like 'title' or 'units'...
            %   variableName = The name of the variable whos attributes you want
            %       to retrieve. If no argument is specified then the
            %       global attributes are returned.
            %
            % Return:
            %   The value associated with the attribute field corresponding to key (and optionally
            %   variableName)
            if nargin < 3
                atlist = obj.attributes;
                val = value4key(atlist, varargin{1});
            elseif nargin == 3
                atlist = obj.attributes(varargin{1});
                val = value4key(atlist, varargin{2});
            else
                warning('NCTOOLBOX:ncugrid:attribute', 'No key or variable specified.');
            end
        end % end attribute

        function gs = getSomeCoordinateData(obj, varName, dataset, meshset)
%             mesh_ident = obj.attribute(varName, 'mesh');
%             mesh_vars = regexp(mesh_ident, ' ', 'split');
%             if nargin < 3
%                 meshset  = obj.meshes.get(str2double(mesh_vars{1}(5:end)-1));
%             end
            connectivity = meshset.getMesh().getTopology().getFaceNodeConnectivityVariable();
            coordinates = obj.attribute(varName, 'coordinates');
            coordinate_vars = regexp(coordinates, ' ', 'split');

             location = obj.attribute(varName, 'location');
%             connectivity_var = obj.attribute([location, '_connectivity'], mesh_vars{1});

            % Get coordinates (time, lat, lon,...)
            for i = 1:length(coordinate_vars)
                switch coordinate_vars{i}
                    case 'time'
                        v = meshset.getMeshVariableByName(coordinate_vars{i});
                        if isempty(v)
                            v = dataset.getNetcdfFile().findVariable(coordinate_vars{i});
                        else
                            %
                        end
                        array = v.read();
                        values = array.copyToNDJavaArray();
                        gs.(coordinate_vars{i}) = values;
                        clear values
                    otherwise
                        v = meshset.getMeshVariableByName(coordinate_vars{i});
                        if isempty(v)
                            v = dataset.getNetcdfFile().findVariable(coordinate_vars{i});
                        else
                            %
                        end
                        array = v.read();
                        values = array.copyToNDJavaArray();
                        gs.(coordinate_vars{i}) = values;
                        clear values
                end
            end % end for

%             switch location
%                 case 'node'
%                     locs = meshset.getMesh().getUniqueNodes();
%                     for i = 1:locs.size();
%                         %                         point = locs.get(i-1).getGeoPoint();
%                         %                 gs.lat(i,1) = point.getLatitude();
%                         %                 gs.lon(i,1) = point.getLongitude();
%                         gs.index(i,1) = locs.get(i-1).getDataIndex();
%                     end
%                 case 'face'
%                     locs = meshset.getMesh().getUniqueFaces();
%                     for i = 1:locs.size();
%                         gs.index(i,1) = locs.get(i-1).getDataIndex();
%                     end
%                     %                 case 'edge'
%                     %                     locs = meshset.getMesh().getUniqueEdges();
%             end
            % Kyle made a method to do this automatically mesh.getNodeLatLons() and getNodeIndexes()
            %             for i = 1:locs.size();
            %                 point = locs.get(i-1).getGeoPoint();
            % %                 gs.lat(i,1) = point.getLatitude();
            % %                 gs.lon(i,1) = point.getLongitude();
            %                 gs.index(i,1) = locs.get(i-1).getDataIndex();
            %             end

            try
                % Get connectivity
                %                     v = obj.netcdfugrid.getNetcdfFile().findVariable(connectivity_var);
%                 v = meshset.getMesh();
%                 v = v.getConnectivityVariable();
%                 v = connectivity.getVariable();

                array = connectivity.read();
                values = array.copyToNDJavaArray();
                gs.connectivity = values';
                switch location
                    case 'node'
                        gs.index = unique(gs.connectivity);
                end
            catch me
                me.throw();
            end

        end % end

        function d = data(obj, variable, first, last, stride)
            % NCUGRID.DATA - Get variable data based on vector of elements (nodes, faces, etc.) representation of the unstructured data.
            %
            % Usage:
            %              d = nc.data(variableName)
            %              d = nc.data(variableName, first)
            %              d = nc.data(variableName, first, last)
            %              d = nc.data(variableName, first, last, stride)
            %

            variable = ucar.nc2.NetcdfFile.escapeName(variable);
            v = obj.netcdfugrid.getNetcdfFile().findVariable(variable);

            if (nargin == 2)
                array = v.read();
                try
                    d = array.copyToNDJavaArray(); % this fails if the variable has no java shape/no dimension was assigned
                catch me1
                    try
                        % TODO (Alex added this code) Where is a file where
                        % this code section gets called?
                        d = array.toString;  % different way to get single value out of java array
                        d = d.toCharArray';  % must transpose
                        d = str2double(d);   % matlab string to matlab numeric
                    catch me2
                        ex = MException('NCTOOLBOX:ncugrid:data', ['Failed to open "' variable '" in ' url]);
                        ex = ex.addCause(me2);
                        ex.throw;
                    end
                end
                d = d;
            else
                s = obj.size(variable);

                % Fill in missing arguments
                % default stride is 1
                if (nargin < 5)
                    stride = ones(1, length(s));
                end

                % Default last is the end
                if (nargin < 4)
                    last = s;
                end

                % Construct the range objects needed to subset the data
                n = max(size(obj.size(variable)));
                ranges = java.util.ArrayList(n);
                for i = 1:n
                    ranges.add(ucar.ma2.Range(first(i) - 1, last(i) - 1, stride(i)));
                end

                array = v.read(ranges);
                d = array.copyToNDJavaArray();
            end
        end % end data

        function ig = grid_interop(src, varName, first, last, stride)
            % NCUGRID.GRID_INTEROP - Method to get the coordinate variables and their data as a
            % a structure with standardized field names for lat, lon, time, and z. Other coordiante variables
            % that are not recognized by netcdf-java as the previous four types have field names directly
            % taken from their variable names.
            % Usage:
            %          >> gridstruct = geovar.grid_interop(1,:,:,1:2:50);
            %
            g = src.grid(varName, first, last, stride);
            names = fieldnames(g);

            for i = 1:length(names); % loop through fields returned by grid
                tempname = names{i};
                javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                try
                    type = char(javaaxisvar.getAxisType);
                    if isempty(type)
                        ig.(tempname) = g.(tempname);
                    else
                        switch type
                            case 'Height'
                                pos_z = char(javaaxisvar.getPositive());
                                if strcmp(pos_z, 'POSITIVE_DOWN')
                                    tmp = g.(tempname);
                                    ig.z = tmp.*-1; %adjust for positive direction
                                else
                                    ig.z = g.(tempname);
                                end

                            case 'GeoZ'
                                pos_z = char(javaaxisvar.getPositive());
                                if strcmp(pos_z, 'POSITIVE_DOWN')
                                    tmp = g.(tempname);
                                    ig.z = tmp.*-1; %adjust for positive direction
                                else
                                    ig.z = g.(tempname);
                                end

                            case 'Time'
                                tmp = g.(tempname);
                                t_converted = src.dataset.time(tempname, tmp);
                                ig.time = t_converted;

                                %                     case 'RunTime'
                                %                       tmp = obj.dataset.data(name, vFirst, vLast, vStride);
                                %                       t_converted = obj.dataset.time(name, tmp);
                                %                       data.(type) = t_converted;

                            case 'Lon'
                                tmp = g.(tempname);
                                ind = find(tmp > 180); % convert 0-360 convention to -180-180
                                tmp(ind) = tmp(ind)-360;
                                ig.lon = tmp;

                            case 'Lat'
                                ig.lat = g.(tempname);

                            otherwise
                                ig.(tempname) = g.(tempname);

                        end % end switch on type
                    end % end is type empty or not if statement
                catch
                    ig.(tempname) = g.(tempname);
                end
            end % end loop through field names

        end % grid_interop end

        function gr = grid(obj, varName, first, last, stride)
            % NCUGRID.GRID - Return a matlab structure of the coordinate variables of 'variableName',
            % subsetted using first, last, and stride index arguments in relation to the variable of interest
            % 'variableName'.
            %
            % Usage:
            %              d = nc.grid(variableName)
            %              d = nc.grid(variableName, first)
            %              d = nc.grid(variableName, first, last)
            %              d = nc.grid(variableName, first, last, stride)
            %
            mesh_ident = obj.attribute(varName, 'mesh');
            mesh_vars = regexp(mesh_ident, ' ', 'split');
            if nargin < 3
                mesh  = obj.meshes.get(str2double(mesh_vars{1}(5:end)-1));
            end

            coordinates = obj.attribute(varName, 'coordinates');
            coordinate_vars = regexp(coordinates, ' ', 'split');

            location = obj.attribute(varName, 'location');
            if strcmp(location, 'node')
                connectivity_var = obj.attribute(mesh_vars{1}, ['face_', location, '_connectivity']);
            else

            end
            % Get coordinates (time, lat, lon,...)
            for i = 1:length(coordinate_vars)
                v = obj.netcdfugrid.getNetcdfFile().findVariable(coordinate_vars{i});
                vs = obj.size(v.getName);
                s = obj.size(varName);
                %         if numel(vs) > 0 % Added to solve work around somedata calls that involve variables with
                %           % no netcdf dim. (This will be frequent in some OOI files.)
                if (length(vs) == length(s))
                    %% case: sizes are the same
                    if isequal(vs, s)
                        vFirst = first;
                        vLast = last;
                        vStride = stride;
                    else
                        me = MException('NCTOOLBOX:ncvariable:somedata', ...
                            ['The data size of the coordinate variable,' ...
                            name ', does not fit the size of ' obj.name]);
                        me.throw;
                    end

                elseif length(vs) == 1
                    %% case: singleton dimension. Find side of data with
                    % the same length

                    % TODO: the following line  will give bogus results if
                    % the data has 2 dimensions of the same length
                    dim = find(s == vs, 1);
                    if ~isempty(dim)
                        vFirst = first(dim);
                        vLast = last(dim);
                        vStride = stride(dim);
                    else
                        me = MException('NCTOOLBOX:ncvariable:somedata', ...
                            ['The data size of the coordinate variable,' ...
                            name ', does not fit the size of ' obj.name]);
                        me.throw;
                    end

                else
                    %% case: variable is coordiantes. Look for size
                    % TODO this is a lame implementation.
                    dim = find(s == vs(1), 1);
                    if ~isempty(dim)
                        for j = 2:length(vs)
                            if vs(j) ~= s(dim + j - 1)
                                me = MException('NCTOOLBOX:ncvariable:somedata', ...
                                    ['The data size of the coordinate variable,' ...
                                    name ', does not fit the size of ' obj.name]);
                                me.throw;
                            end
                        end
                        k = dim:dim + length(vs) - 1;
                        vFirst = first(k);
                        vLast = last(k);
                        vStride = stride(k);
                    end
                end

                n = max(size(obj.size(v.getName)));
                ranges = java.util.ArrayList(n);
                for j = 1:n
                    ranges.add(ucar.ma2.Range(vFirst(j) - 1, vLast(j) - 1, vStride(j)));
                end

                array = v.read(ranges);
                values = array.copyToNDJavaArray();
                gr.(coordinate_vars{i}) = values;
                clear values

            end
            %         else
            % not needed i dont think...

            % Get connectivity
            v = obj.netcdfugrid.getNetcdfFile().findVariable(connectivity_var);
            array = v.read();
            values = array.copyToNDJavaArray();
            gr.connectivity = values;

        end % end grid

        function vs = size(obj, varName)
            % NCUGRID.SIZE - Return the dimensional sizes of the variable of interest.
            %
            % Usage:
            %              d = nc.size(variableName)
            %
            v = obj.dataset.netcdf.findVariable(varName);
            vs = v.getShape;
            vs = vs';
        end

        function d = timewindowij(src, varargin)
            % NCGEOVARIABLE.TIMEWINDOWIJ - Function to get indices from start and stop times for sub-
            % setting. TODO: There must be a better/fast way to do this using the java library.
            % Useage: >> timestruct = nc.timewindowij([2004 1 1 0 0 0], [2005 12 31 0 0 0]);
            %              >> timestruct = nc.timewindowij(731947, 732677);
            % Result: time.time = time data
            %            time.index = time indices
            s = src.size(varargin{1});
            first = ones(1, length(s));
            last = s;
            stride = first;
            g = src.grid_interop(first, last, stride);

            if isfield(g, 'time') % are any of the fields recognized as time explictly
                starttime = datenum(varargin{2});
                stoptime = datenum(varargin{3});
                if isempty(starttime)
                    starttime = g.time(1);
                end
                if isempty(stoptime)
                    stoptime = g.time(end);
                end

                t_index1 = g.time >= starttime;
                t_index2 = g.time <= stoptime;
                d.index = find(t_index1==t_index2);
                d.time = g.time(d.index);
            else
                me = MException(['NCTOOLBOX:ncgeovariable:timewindowij'], ...
                    'No grid variable returned as time.');
                me.throw;
            end
        end % end timewindowij

    end % end methods

end % end class

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function test_example_fem

% WALLTIME 08:00:00
% MEM 12gb
% DEPENDENCY ft_prepare_headmodel ft_prepare_mesh ft_datatype_segmentation

[ftver, ftpath] = ft_version;

%% read mri
mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri'));

%% segmentation
cfg          = [];
cfg.output   = {'gray', 'white', 'csf', 'skull', 'scalp'};
segmentedmri = ft_volumesegment(cfg, mri);

%% mesh
cfg        = [];
cfg.shift  = 0.3;
cfg.method = 'hexahedral';
mesh       = ft_prepare_mesh(cfg,segmentedmri);

%% volume conductor
cfg              = [];
cfg.method       = 'simbio';
cfg.conductivity = [0.33 0.14 1.79 0.01 0.43];
headmodel        = ft_prepare_headmodel(cfg, mesh);

%% electrode alignment
elec = ft_read_sens(fullfile(ftpath,'template/electrode/standard_1020.elc'));

nas = mri.hdr.fiducial.head.nas;
lpa = mri.hdr.fiducial.head.lpa;
rpa = mri.hdr.fiducial.head.rpa;

%
fiducials.pos     = [nas; lpa; rpa];
fiducials.label   = {'Nz','LPA','RPA'};
fiducials.unit    = 'mm';

cfg          = [];
cfg.method   = 'fiducial';
cfg.target   = fiducials;
cfg.elec     = elec;
cfg.fiducial = {'Nz', 'LPA', 'RPA'};
elec_align   = ft_electroderealign(cfg);

% add 12 mm to x-axis
n=size(elec_align.chanpos,1);
for i=1:n
 elec_align.chanpos(i,1)=elec_align.chanpos(i,1)+12;
 elec_align.elecpos(i,1)=elec_align.elecpos(i,1)+12;
end

%% make the sourcemodel/grid

% At the moment the sourcemodel is defined prior
% to the leadfield because ft_prepare_sourcemodel does not automatically create
% a sourcemodel based on a hexahedral headmodel.

cfg                 = [];
cfg.mri             = mri;
cfg.resolution      = 3;
sourcemodel         = ft_prepare_sourcemodel(cfg);
sourcemodel         = ft_convert_units(sourcemodel, headmodel.unit);
sourcemodel.inside(500:end) = false; % do only a few dipoles

cfg            = [];
cfg.headmodel  = headmodel;
cfg.elec       = elec_align;
cfg.sourcemodel = sourcemodel;
cfg.channel    = elec_align.label(4:6); % do only a few channels, to save time
lf             = ft_prepare_leadfield(cfg);

m = cellfun(@mean,lf.leadfield(lf.inside),'uniformoutput',false)';
m = cat(1, m{:});
assert(all(m(:)<eps));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
 function output = ir_radon_zwart_powell(theta, rr)
%function output = ir_radon_zwart_powell(theta, yr)
%|
%| Compute analytic 2D Radon transform of Zwart-Powell box spline.
%|
%| in
%| theta	ray angle in radian
%| rr	distance between the point and the ray (normalized by the pixel size)
%| output: radon transform of the Zwart-Powell box spline element
%|
%| This is the modified version of the code written by [1]
%| to avoid symbolic math operation.
%| Code written by Seongjin Yoon, Univ. of Michigan, Jan 2015
%|
%| Reference
%| [1] A. Entezari, M. Nilchian, and M. Unser, "A box spline calculus
%| for the discretization of computed tomography reconstruction problems,"
%| IEEE Trans. Med. Imaging, vol. 31, no. 8, pp. 1532–1541, 2012.
%|
%| 2015-08-10 Jeff Fessler, added self test and parallelized

if nargin == 1 && streq(theta, 'test'), ir_radon_zwart_powell_test, return, end
if nargin < 2, ir_usage, end

dim = size(theta);
theta = theta(:);

zeta = zeros(numel(theta), 4);
zeta(:,1) = cos(theta);
zeta(:,2) = sin(theta);
zeta(:,3) = zeta(:,1) + zeta(:,2);
zeta(:,4) = zeta(:,1) - zeta(:,2);

cond = abs(zeta) >= eps('single');
N = sum(cond,2);

output = BoxSp4(rr(:), zeta, cond, N) ./ factorial(N-1);
output = reshape(output, dim);

function output = BoxSp0(y, N)
%output = heaviside(y) .* y.^(N-1);
if any(size(y) ~= size(N)), keyboard, end
output = (y >= 0) .* y.^(N-1);

function output = BoxSp1(y, zeta, cond, N)
good = cond(:,1);
output = (BoxSp0(y+0.5*zeta(:,1), N) ...
	- BoxSp0(y-0.5*zeta(:,1), N)) ./ zeta(:,1);
output(~good) = BoxSp0(y(~good), N(~good));

function output = BoxSp2(y, zeta, cond, N)
good = cond(:,2);
output = (BoxSp1(y+0.5*zeta(:,2), zeta, cond, N) ...
	- BoxSp1(y-0.5*zeta(:,2), zeta, cond, N)) ./ zeta(:,2);
output(~good) = BoxSp1(y(~good), zeta(~good,:), cond(~good,:), N(~good));

function output = BoxSp3(y, zeta, cond, N)
good = cond(:,3);
output = (BoxSp2(y+0.5*zeta(:,3), zeta, cond, N) ...
	- BoxSp2(y-0.5*zeta(:,3), zeta, cond, N)) ./ zeta(:,3);
output(~good) = BoxSp2(y(~good), zeta(~good,:), cond(~good,:), N(~good));

function output = BoxSp4(y, zeta, cond, N)
good = cond(:,4);
output = (BoxSp3(y+0.5*zeta(:,4), zeta, cond, N) ...
	- BoxSp3(y-0.5*zeta(:,4), zeta, cond, N)) ./ zeta(:,4);
output(~good) = BoxSp3(y(~good), zeta(~good,:), cond(~good,:), N(~good));

% ir_radon_zwart_powell_test()
function ir_radon_zwart_powell_test
theta = linspace(0,pi,181);
r = linspace(-1,1,101) * 2;
[tt rr] = ndgrid(theta, r);
sino = ir_radon_zwart_powell(tt, rr);
im('colorneg', r, theta, sino'), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% ZEROAXES Redraws the axes at the zero values.
% See also: CENTAXES
% Author: Andrew Knight

% Modified February 1994 to omit the zero label if necessary.
% Modified May 12th 2008 to allow  Y axis invertion, by Javier Lopez-Calderon (JLC)
%
% 03-31-2009 : Bug1: incorrect position of x-axis when y-axis is inverted (abs(ymin) ~= abs(ymax)).
% Thanks to Marcus Heldmann.
%
% 04-01-2009 : Bug1 fixed. JLC

function zeroaxes

holdwason = ishold;
axesin    = get(gcf,'CurrentAxes');
bdownf    = get(gcf,'ButtonDownFcn'); % JLC, May 12th 2008
pos       = get(axesin,'position');
ax        = axis;
xscale    = get(axesin,'XScale');
yscale    = get(axesin,'YScale');
xxticks   = get(axesin,'XTick');       %JLC, March 10, 2011
yyticks   = get(axesin,'YTick');       %JLC, March 10, 2011
xmt       = get(axesin,'XMinorTick');  %JLC, March 10, 2011
ymt       = get(axesin,'YMinorTick');  %JLC, March 10, 2011

sentido   = get(axesin,'YDir');       % JLC, May 12th 2008
set(axesin,'visible','off');
xmin = ax(1);
xmax = ax(2);
ymin = ax(3);
ymax = ax(4);
ticklength  = get(axesin,'TickLength');
XAxisHeight = ticklength(1);
YAxisWidth  = ticklength(1);
f   = polyfit([ax(1) ax(2)],[pos(1) pos(1)+pos(3)],1);
XAxisXLimits = polyval(f,[xmin xmax]);
YAxisXLimits = polyval(f,[0 YAxisWidth*abs(xmax - xmin)]);
f   = polyfit([ax(3) ax(4)],[pos(2) pos(2)+pos(4)],1);
YAxisYLimits = polyval(f,[ymin ymax]);
XAxisYLimits = polyval(f,[0 XAxisHeight*abs(ymax - ymin)]);

%
% right position of the x-axis in case of inverted y-axis
% 04-01-2009  JLC
if strcmp(sentido, 'reverse')
      Xaxis_y = 2*pos(2)+pos(4)-XAxisYLimits(1);
else
      Xaxis_y = XAxisYLimits(1);
end

XAxisPosition = [XAxisXLimits(1)
      Xaxis_y
      XAxisXLimits(2) - XAxisXLimits(1)
      XAxisYLimits(2) - XAxisYLimits(1)];

bgcolour = get(gcf,'color');

hX = axes('position',XAxisPosition,...
      'XLim',[xmin xmax],...
      'box','off',...
      'YTick',[],...
      'TickDir','out',...
      'XScale',xscale,...
      'YColor',bgcolour,...
      'color','none');

YAxisPosition = [YAxisXLimits(1)
      YAxisYLimits(1)
      YAxisXLimits(2) - YAxisXLimits(1)
      YAxisYLimits(2) - YAxisYLimits(1)];

hY = axes('position',YAxisPosition,...
      'YLim',[ymin ymax],...
      'box','off',...
      'Xtick',[],...
      'TickDir','out',...
      'YScale',yscale,...
      'XColor',bgcolour,...
      'color','none',...
      'YDir', sentido); % JLC, May 12th 2008

%
% Ticks  JLC, March 10, 2011
%
set(hX,'XTick',xxticks)
set(hY,'YTick',yyticks)
set(hX,'XMinorTick',xmt)
set(hY,'YMinorTick',ymt)

% Get rid of the zero ticks if necessary:
if ymin<0 && ~strcmp(xscale,'log')
      xticks = get(hX,'XTick');
      xticks(xticks==0) = [];
      set(hX,'XTick',xticks)
end

if xmin<0 && ~strcmp(yscale,'log')
      yticks = get(hY,'YTick');
      yticks(yticks==0) = [];
      set(hY,'YTick',yticks)
end

set(gcf,'CurrentAxes',axesin)
set(gcf,'ButtonDownFcn',bdownf)

if ~holdwason
      set(hX,'NextPlot','Replace')
      set(hY,'NextPlot','Replace')
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function togglegrid(dffig)
%TOGGLEGGRID Toggle x and y axes grid on or off

%   $Revision: 1.1.6.2 $  $Date: 2004/01/24 09:35:54 $
%   Copyright 2003-2004 The MathWorks, Inc.

% Get new state -- note uimenu state reflects old state, and
% uitoggletool state reflects new state
h = gcbo;
if isequal(get(h,'Type'),'uimenu')
   onoff = on2off(get(h,'Checked'));
else
   onoff = get(h,'State');
end
dfgetset('showgrid',onoff);

% Change grid
ax = findall(dffig,'Type','axes');
for j=1:length(ax)
   if ~isequal(get(ax(j),'Tag'),'legend')
      set(ax(j),'xgrid',onoff,'ygrid',onoff)
   end
end

% Change menu state
h = findall(dffig,'Type','uimenu','Tag','showgrid');
if ~isempty(h), set(h,'Checked',onoff); end

% Change button state
h = findall(dffig,'Type','uitoggletool','Tag','showgrid');
if ~isempty(h), set(h,'State',onoff); end

author: zouchuhang

repo: LayoutNet


text:
function fiberPath = dtiFiberTensorlinePreinterp(seedPoint, dt6, faImg, voxSize, faThresh, angleThresh, ...
    stepSizeMm, wPuncture, whichAlgorithm)

% fiberPath = dtiFiberTensorlinePreinterp(seedPoint, dt6, faImg, voxSize, faThresh, angleThresh, stepSizeMm, wPuncture, whichAlgorithm)
%
% Implementation of Lazar's Tensorline algorithm. Eg:
% White Matter Tractography Using Diffusion Tensor Deflection.
% Lazar, et. al. Human Brain Mapping 18:306 321(2003).
%
%
% Can be used to find Tensorline, FACT, or TEND fiber tracing
% Choosing Algorithm using whichAlgorithm:
%   0 = TensorLine
%   1 = FACT
%   2 = TEND
%
% Step Size:
% Determines resolution of DTI grid; default is 2x2x2mm voxels
% 1. Step Size < voxel size of inputted DTI matrix.
% Algorithm is forced to reinterpolate tensors at each point (SLOW!).
% 2. Step Size > voxel size of inputted matrix
% Step Size is reinitialized to voxel size of inputted matrix.
% Algorithm uses "nearest neighbor" concept to find proper tensor
% from the inputted dt6 file
%
% RETURNS:
% fiberPath, a list of coordinates (in voxels, but may be real-valued)
% that define the fiber path.
%
% HISTORY:
%   2004.02.09 GSM (gmulye@stanford.edu) wrote it.
%   2004.02.09 GSM Minor updates, more userdefined parameters
%   2004.02.13 GSM Modified to use pre-interpolated tensors

%voxSize = [1 1 1]'; %%%FORCE SMALL VOXELS
%oneMmDt6GridSize = [145 196 131 6]; %Size of dt6 array with 1x1x1 voxels
%origDt6GridSize = size(dt6); %Size of dt6 array inputted into function
%origDt6VoxSize = (origDt6GridSize + 1) ./ oneMmDt6GridSize; %vox size of inputted array

stepSize = stepSizeMm*ones(3,1);

if all(voxSize <= stepSize)
    reinterp = 0;
else
    reinterp = 1;
end

%stepSize = stepSizeMm./voxSize; %Step size in voxels - voxSize is 3x1 column vector

% Initialize variables for tracing
% Find major direction (originalDir) of seedPoint voxel
if (reinterp == 0)
    %   voxCoords = nearestNeighbor(seedPoint,voxSize);
    voxCoords = round(seedPoint);
    dt6Tensor = dt6(voxCoords(1),voxCoords(2),voxCoords(3),1:6);
    dt6Tensor = dt6Tensor(:); %Vectorize
elseif (reinterp == 1)
    dt6Tensor = findTensor(seedPoint,dt6,voxSize);
end
[originalDir junk] = majorEigVec(dt6Tensor);
nextDir = originalDir; % Initialize to seedPoint direction
nextPosition = seedPoint; %Absolute position in voxels
fiberPath = seedPoint; % First point is seed point
fiberPath = tracer(whichAlgorithm,nextPosition,originalDir,voxSize,dt6,reinterp,faImg,faThresh,angleThresh,wPuncture,stepSize,fiberPath,1);
fiberPath = tracer(whichAlgorithm,nextPosition,originalDir,voxSize,dt6,reinterp,faImg,faThresh,angleThresh,wPuncture,stepSize,fiberPath,-1);
return;

function fiberPath = tracer(whichAlgorithm,nextPosition,nextDir,voxSize,dt6,reinterp,faImg,faThresh,angleThresh,wPuncture,stepSize,fiberPath,fwdBkwd)
%==================
% MAIN FUNCTION
%==================
%Traces fiberpath forwards and backwards
iter = 0;
done = 0;
maxIter = 1000;
imSize = size(dt6(:,:,:,1));
nextDir = fwdBkwd*nextDir;
while (~done & iter<maxIter)
    dir = nextDir; %Direction vectors are all 3x1 column vectors
    currentPosition = nextPosition; %Position coordinates are 1x3 row vectors
    if(any(currentPosition<1) | any(currentPosition>imSize))
        disp('Tracking terminated: path wandered outside image data.');
        done = 1;
    else
        % Get the FA for this voxel
        vCoord = floor(currentPosition);
        fa = faImg(vCoord(1), vCoord(2), vCoord(3));

        % check dir data is valid (we are in range)
        if (fa <= faThresh | sum(dir) == 0 | isnan(fa))
            disp(['Tracking terminated: fa=',num2str(fa)]);
            done = 1;
        else
            nextPosition = currentPosition + (stepSize'./voxSize') .*dir'; %convert stepSize to voxels
            % Find next diffusion tensor
            if (reinterp == 0)
                %voxCoords = nearestNeighbor(nextPosition,voxSize);
                voxCoords = round(nextPosition);
                if(any(voxCoords==0) | any(voxCoords>imSize) | ~any(isfinite(voxCoords)))
                    disp(['Tracking terminated: path wandered outside image']);
                    done = 1;
                    voxCoords(voxCoords==0) = 1;
                else
                    nextDt6Tensor = dt6(voxCoords(1),voxCoords(2),voxCoords(3),1:6);
                    nextDt6Tensor = nextDt6Tensor(:); %Vectorize
                end
            elseif (reinterp == 1)
                nextDt6Tensor = findTensor(nextPosition,dt6,voxSize);
            end
            % FACT direction, checked to see if angle is less than
            % threshold
            [factDir, nextTensor] = majorEigVec(nextDt6Tensor);
            [smallAngle,direction,angle] = checkAngle(dir,factDir,angleThresh);
            if (smallAngle & (direction == -1))
                factDir = -factDir;
            end
            % TEND direction
            tendDir = nextTensor*dir;
            nd = norm(tendDir);
            if(nd<=0)
                tendDir = tendDir/nd;
            else
                tendDir = NaN;
            end
            % Tensorline Algorithm
            if (whichAlgorithm == 0)
                if(tendDir==NaN)
                    disp(['Tracking terminated: tendDir = NaN']);
                    done = 1;
                else
                    if (~smallAngle) % Turning too sharp, ignore FACT vector
                        nextDir = (1-wPuncture)*dir + wPuncture*tendDir;
                    else  % Normal case: Use FACT vector
                        nextDir = fa*factDir + (1-fa)*((1-wPuncture)*dir + wPuncture*tendDir);
                    end
                end
                % FACT only
            elseif (whichAlgorithm == 1)
                if (~smallAngle) % Angle too big
                    done = 1; % End tracing due to sharp turn
                    disp(['Tracking terminated: angle (' num2str(round(angle)) ') exceeds threshold.']);
                else
                    nextDir = factDir;
                end
                % TEND only
            elseif (whichAlgorithm == 2)
                if(tendDir==NaN)
                    disp(['Tracking terminated: tendDir = NaN']);
                    done = 1;
                else
                    nextDir = tendDir;
                end
            end
            % Append to fiberpath
            if (fwdBkwd == 1)
                fiberPath = [fiberPath; nextPosition];
            elseif (fwdBkwd == -1)
                fiberPath = [nextPosition;fiberPath];
            end
        end
    end
end
return;

function dt6Tensor = findTensor(point,dt6,voxSize)
persistent initialized;
% Finds tensor interpolation at any point and return dt6 tensor
curPosMm = point.*voxSize';
if(isempty(initialized))
    dt6Tensor = dtiTensorInterp(dt6, curPosMm, voxSize', 1);
    initialized = 1;
else
    dt6Tensor = dtiTensorInterp([], curPosMm, voxSize', 1);
end
return;

function intVox = nearestNeighbor(realCoord,voxSize)
% Rounds mm input into voxel coordinates
decVox = realCoord./voxSize'; %decimal voxels is 1x3
intVox = round(decVox);
return;

function [majorDir, fullTensor] = majorEigVec(dt6Tensor)
% Finds major direction for given tensor (in 3x3 format)
fullTensor = [dt6Tensor(1) dt6Tensor(4) dt6Tensor(5); ...
        dt6Tensor(4) dt6Tensor(2) dt6Tensor(6); ...
        dt6Tensor(5) dt6Tensor(6) dt6Tensor(3)];
[eigVec,eigVal] = eig(fullTensor);
[maxVal,i] = max(max(eigVal));
majorDir = eigVec(:,i);
%majorDir = [majorDir(2) majorDir(3) majorDir(1)]';
return;

function [angleCheck, direction, angle] = checkAngle(a,b,angleThresh)
% Checks the angle between 2 vectors, returns if angle is less than thresh
% The angle between two vectors is given by acos(aDOTb/{mag(a)*mag(b)})
anglePos = 180*acos(a'*b)/pi; % In degrees; both vectors are unit vectors
angleNeg = 180*acos(-a'*b)/pi; % Angle with one vector reversed
[angle,i] = min(abs([anglePos angleNeg]));
direction = 0;
if (angle <= angleThresh)
    angleCheck = 1; % Angle is within permissable range
    if (i==2)
        direction = -1; % Use FACT vector in opposite direction
    end
else
    angleCheck = 0; % Does not pass check
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [concentrationMatrix, excRxnNames, timeVec, biomassVec] = dynamicFBA(model, substrateRxns, initConcentrations, initBiomass, timeStep, nSteps, plotRxns, exclUptakeRxns)
% Performs dynamic FBA simulation using the static optimization approach
%
% USAGE:
%
%    [concentrationMatrix, excRxnNames, timeVec, biomassVec] = dynamicFBA(model, substrateRxns, initConcentrations, initBiomass, timeStep, nSteps, plotRxns, exclUptakeRxns)
%
% INPUTS:
%    model:                  COBRA model structure
%    substrateRxns:          List of exchange reaction names for substrates
%                            initially in the media that may change (e.g. not
%                            h2o or co2)
%    initConcentrations:     Initial concentrations of substrates (in the same
%                            structure as `substrateRxns`)
%    initBiomass:            Initial biomass (must be non zero)
%    timeStep:               Time step size
%    nSteps:                 Maximum number of time steps
%
% OPTIONAL INPUTS:
%    plotRxns:               Reactions to be plotted (Default = {'EX_glc(e)', 'EX_ac(e)', 'EX_for(e)'})
%    exclUptakeRxns:         List of uptake reactions whose substrate concentrations do not change
%                            (Default = {'EX_co2(e)', 'EX_o2(e)', 'EX_h2o(e)', 'EX_h(e)'})
%
% OUTPUTS:
%    concentrationMatrix:    Matrix of extracellular metabolite concentrations
%    excRxnNames:            Names of exchange reactions for the EC metabolites
%    timeVec:                Vector of time points
%    biomassVec:             Vector of biomass values
%
% If no initial concentration is given for a substrate that has an open
% uptake in the model (i.e. `model.lb < 0`) the concentration is assumed to
% be high enough to not be limiting. If the uptake rate for a nutrient is
% calculated to exceed the maximum uptake rate for that nutrient specified
% in the model and the max uptake rate specified is > 0, the maximum uptake
% rate specified in the model is used instead of the calculated uptake
% rate.
%
% NOTE:
%
%    The dynamic FBA method implemented in this function is essentially
%    the same as the method described in
%    [`Varma, A., and B. O. Palsson. Appl. Environ. Microbiol. 60:3724 (1994)`].
%    This function does not implement the dynamic FBA using dynamic optimization approach
%    described in [`Mahadevan, R. et al. Biophys J, 83:1331-1340 (2003)`].
%
% .. Author: - Markus Herrgard 8/22/06

global WAITBAR_TYPE

if (nargin < 7)
    plotRxns = {'EX_glc(e)','EX_ac(e)','EX_for(e)'};
end

% Uptake reactions whose substrate concentrations do not change
if (nargin < 8)
    exclUptakeRxns = {'EX_co2(e)','EX_o2(e)','EX_h2o(e)','EX_h(e)'};
end

% Find exchange rxns
excInd = findExcRxns(model,false);
excInd = excInd & ~ismember(model.rxns,exclUptakeRxns);
excRxnNames = model.rxns(excInd);
length(excRxnNames)
% Figure out if substrate reactions are correct
missingInd = find(~ismember(substrateRxns,excRxnNames));
if (~isempty(missingInd))
    for i = 1:length(missingInd)
        fprintf('%s\n',substrateRxns{missingInd(i)});
    end
    error('Invalid substrate uptake reaction!');
end

% Initialize concentrations
[~, substrateMatchInd] = ismember(substrateRxns,excRxnNames);
concentrations = zeros(length(excRxnNames),1);
concentrations(substrateMatchInd) = initConcentrations;

% Deal with reactions for which there are no initial concentrations
originalBound = -model.lb(excInd);
noInitConcentration = (concentrations == 0 & originalBound > 0);
concentrations(noInitConcentration) = 1000;

biomass = initBiomass;

% Initialize bounds
uptakeBound =  concentrations/(biomass*timeStep);

% Make sure bounds are not higher than what are specified in the model
aboveOriginal = (uptakeBound > originalBound) & (originalBound > 0);
uptakeBound(aboveOriginal) = originalBound(aboveOriginal);
model.lb(excInd) = -uptakeBound;

concentrationMatrix = sparse(concentrations);
biomassVec = biomass;
timeVec(1) = 0;

fprintf('Step number\tBiomass\n');
showprogress(0,'Dynamic FBA analysis in progress ...');
for stepNo = 1:nSteps
    % Run FBA
    sol = optimizeCbModel(model,'max','one');
    mu = sol.f;
    if (sol.stat ~= 1 || mu == 0)
        fprintf('\nNo feasible solution - nutrients exhausted. Biomass:\t %f\n', biomass);
        break;
    end
    uptakeFlux = sol.x(excInd);
    biomass = biomass*exp(mu*timeStep);
    %biomass = biomass*(1+mu*timeStep);
    biomassVec(end+1) = biomass;

    % Update concentrations
    concentrations = concentrations - uptakeFlux/mu*biomass*(1-exp(mu*timeStep));
    %concentrations = concentrations + uptakeFlux*biomass*timeStep;
    concentrations(concentrations <= 0) = 0;
    concentrationMatrix(:,end+1) = sparse(concentrations);

    % Update bounds for uptake reactions
    uptakeBound =  concentrations/(biomass*timeStep);
    % This is to avoid any numerical issues
    uptakeBound(uptakeBound > 1000) = 1000;
    % Figure out if the computed bounds were above the original bounds
    aboveOriginal = (uptakeBound > originalBound) & (originalBound > 0);
    % Revert to original bounds if the rate was too high
    uptakeBound(aboveOriginal) = originalBound(aboveOriginal);
    uptakeBound(abs(uptakeBound) < 1e-9) = 0;

    model.lb(excInd) = -uptakeBound;

    if WAITBAR_TYPE ~= 1
        fprintf('%d\t%f\n',stepNo,biomass);
    end
    showprogress(stepNo/nSteps);
    timeVec(stepNo+1) = stepNo*timeStep;
end

selNonZero = any(concentrationMatrix>0,2);
concentrationMatrix = concentrationMatrix(selNonZero,:);
excRxnNames = excRxnNames(selNonZero);
selPlot = ismember(excRxnNames,plotRxns);

% Plot concentrations as a function of time
clf
subplot(1,2,1);
plot(timeVec,biomassVec);
axis tight
title('Biomass');
subplot(1,2,2);
plot(timeVec,concentrationMatrix(selPlot,:));
axis tight
legend(strrep(excRxnNames(selPlot),'EX_',''));

author: opencobra

repo: cobratoolbox


text:
function sftpack_test04 ( )

%*****************************************************************************80
%
%% TEST04 tests R8VEC_SQCTB and R8VEC_SQCTF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
  n = 256;

  ahi = 5.0;
  alo = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  For real slow quarter wave cosine transforms,\n' );
  fprintf ( 1, '  R8VEC_SQCTF does a forward transform;\n' );
  fprintf ( 1, '  R8VEC_SQCTB does a backward transform.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of data items is N = %d\n', n );
%
%  Set the data values.
%
  seed = 123456789;

  [ x, seed ] = r8vec_uniform ( n, alo, ahi, seed );

  r8vec_print_part ( n, x, 10, '  The original data:' );
%
%  Compute the coefficients.
%
  y = r8vec_sqctf ( n, x );

  r8vec_print_part ( n, y, 10, '  The cosine coefficients:' );
%
%  Now compute inverse transform of coefficients.  Should get back the
%  original data.

  x = r8vec_sqctb ( n, y );

  r8vec_print_part ( n, x, 10, '  The retrieved data:' );

  return
end


text:
function theta = angle2Points(varargin)
%ANGLE2POINTS Compute horizontal angle between 2 points
%
%   ALPHA = angle2Points(P1, P2),
%   Pi are either [1*2] arrays, or [N*2] arrays, in this case ALPHA is a
%   [N*1] array. The angle computed is the horizontal angle of the line
%   (P1 P2)
%   Result is always given in radians, between 0 and 2*pi.
%
%   See Also:
%   points2d, angles2d, angle3points, normalizeAngle, vectorAngle
%
%
% ---------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% created the 02/03/2007.
% Copyright 2010 INRA - Cepia Software Platform.

%   HISTORY:
%   2011-01-11 use bsxfun

% process input arguments
if length(varargin)==2
    p1 = varargin{1};
    p2 = varargin{2};
elseif length(varargin)==1
    var = varargin{1};
    p1 = var(1,:);
    p2 = var(2,:);
end

% ensure data have correct size
n1 = size(p1, 1);
n2 = size(p2, 1);
if n1~=n2 && min(n1, n2)>1
    error('angle2Points: wrong size for inputs');
end

% angle of line (P2 P1), between 0 and 2*pi.
dp = bsxfun(@minus, p2, p1);
theta = mod(atan2(dp(:,2), dp(:,1)) + 2*pi, 2*pi);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%Normalizes homogenous coordinates such that the last coordinate is 1
%You can use any dimension of the vectors
%
%Author:    Christian Wengert,
%           Institute of Computer Vision
%           Swiss Federale Institute of Technology, Zurich (ETHZ)
%           wengert@vision.ee.ethz.ch
%           www.vision.ee.ethz.ch/~cwengert/
%
%Input:     x       unnormalized homogenous coordinates
%
%Output     y       normalized homogenous coordinates
%
%Syntax:    y = normalizeHomogenousCoordinates(x)

function y = normalizeHomogenousCoordinates(x)

    %get dimension of array
    ni = size(x,1);
    nj = size(x,2);
    %go through
    for j=1:1:nj
        y(:,j) = x(:,j)./x(ni,j);
    end
    y(ni,:) = ones(1,nj);

author: christianwengert


text:
classdef CurveBall < solvers.Solver
  % CurveBall - Newton solver variant
  %
  % This variant works by 2 gradient descent updates: one for the linear
  % system (inversion), and another for the main objective.
  %
  % Assumptions:
  % - Network's output is a single tensor of predictions (not loss).
  % - Forward pass was executed just before calling this solver's step().

  properties
    loss = 'logistic'  % loss to use: 'ls' (least-squares), 'logistic' (softmax-log)

    lambda = 1  % Hessian regularization (0 = pure Newton)

    momentum = 0.9  % forget factor for linear system solution
    beta = 0.1  % learning rate for linear system GD step

    autoparam = true  % adapt momentum and beta
    autoparam_interval = 1  % only run autoparam every N steps. NOTE: autolambda depends on autoparam running in the same step.
    autoparam_print = false  % to print momentum/beta for debugging

    autolambda = true  % adapt lambda
    autolambda_w1 = 0.999  % factor, how much to change lambda each time
    autolambda_interval = 5  % only run autolambda every N steps. NOTE: autolambda depends on autoparam running in the same step.
    autolambda_thresh = 0.5  % rho thresholds to trigger lambda change; [low, high] vector, or scalar T to use [T, 2 - T]

    labels  % labels for this mini-batch (batch size assumed to be last dim.; cannot run with batch of 1)
    net  % pointer to network object

    count = 0  % iteration count

    state  % store for momentum (z)

    ratio = 0  % autolambda ratio is stored here for debugging
    last_loss  % last computed loss value
  end

  methods
    function o = CurveBall(varargin)
      % parse generic Solver arguments
      % default weight decay is 0, but can override with name-value pair
      varargin = o.parseGenericArgs([{'weightDecay', 0}, varargin]);

      % parse arguments specific to this solver
      vl_parseprop(o, varargin, {'loss', 'lambda', 'beta', 'momentum', ...
        'autoparam', 'autoparam_interval', 'autoparam_print', 'autolambda', ...
        'autolambda_w1', 'autolambda_thresh', 'autolambda_interval'});
      % don't let parent class clear some parameters
      o.conserveMemory = false;
      o.state = {};
    end

    function w = gradientStep(o, w, ~, lr, decay)
      % store properties in local vars
      net = o.net;  %#ok<*PROPLC>
      lambda = o.lambda;

      % initialize state if needed
      z = o.state;
      if isempty(z)
        z = cell(size(w));
        for i = 1:numel(w)
          z{i} = zeros(size(w{i}), 'like', w{i});
        end
      end

      use_decay = any(decay);
      assert(~use_decay, 'Not implemented.');

      %
      % do two GD steps. first update z (linear system state), then use the
      % result (whitened gradient estimate) to update the parameters w.
      %
      % zdelta = J' * Hl * J * z + lambda * z + J' * Jl'
      %
      % znew = momentum * z - beta * zdelta
      %
      % wnew = w + lr * znew
      %

      % compute factor Jz = J * z (for Hessian term)
      Jz = fmad(net, z);

      % evaluate loss, loss gradient Jl, and loss Hessian (Hl).
      % compute hessian_term = J' * Hl * Jz + Jl
      [Jz_, Jl, o.last_loss] = hessian_grad_loss(net, o.loss, o.labels, Jz);

      % backpropagate Jz_ + Jl
      hessian_term = backward(net, Jz_ + Jl);

      % compute z update (gradient descent on linear system)
      delta_z = cell(size(w));
      for i = 1:numel(w)
        delta_z{i} = hessian_term{i} + lambda * z{i};
      end

      %
      % automatic hyperparameters rho (momentum) and beta
      %

      if o.autoparam && rem(o.count, o.autoparam_interval) == 0
        % compute factor Jz = J * z (for Hessian term)
        Jdz = fmad(net, delta_z);

        % evaluate loss, loss gradient Jl, and loss Hessian (Hl).
        % compute hessian_term = J' * Hl * Jz
        Jdz_ = hessian_grad_loss(net, o.loss, o.labels, Jdz);

        % compute momentum and update by solving a 2x2 system A * x = b
        A11 = gather(Jdz(:)' * Jdz_(:));
        A12 = gather(Jz(:)' * Jdz_(:));
        A22 = gather(Jz(:)' * Jz_(:));

        b1 = gather(Jl(:)' * Jdz(:));
        b2 = gather(Jl(:)' * Jz(:));

        for i = 1:numel(w)
          % compute the system we want to invert
          z_vec = z{i}(:);
          dz_vec = delta_z{i}(:);

          % expand scalar if needed
          if isscalar(z_vec), z_vec(1:numel(dz_vec),1) = z_vec; end

          A11 = A11 + gather(dz_vec' * dz_vec) * lambda;  % accumulating scalars is faster (JIT)
          A12 = A12 + gather(dz_vec' * z_vec) * lambda;
          A22 = A22 + gather(z_vec' * z_vec) * lambda;
        end

        % compute beta and momentum coefficient
        A = double([A11, A12; A12, A22]);
        b = double([b1; b2]);
        m_b = A \ b;
        beta = m_b(1);
        momentum = -m_b(2);

        % sanity check
        if ~isfinite(momentum), momentum = 0; end
        if ~isfinite(beta), beta = 0; end

        % store values for next iterations if needed
        o.momentum = momentum;
        o.beta = beta;
      else
        momentum = o.momentum;
        beta = o.beta;
      end

      if o.autoparam_print
        fprintf('mom: %.1g beta: %.1g ', momentum, beta);
      end

      % update parameters with computed beta and momentum
      for i = 1:numel(w)
        % update linear system state
        z{i} = vl_taccum(momentum, z{i}, -beta, delta_z{i}) ;

        % update parameters (returned as this method's output)
        w{i} = vl_taccum(1, w{i}, lr(i), z{i}) ;
      end
      o.state = z;

      %
      % automatic lambda (trust region) hyperparameter
      %

      if o.autolambda && rem(o.count, o.autolambda_interval) == 0
        h_curr = o.last_loss ;
        M = -0.5 * m_b(:)' * b(:);

        % compute new value
        labels = o.labels;
        if isvector(labels)  % categorical labels
          batch_size = numel(labels);
        else  % regression targets, get last dimension
          batch_size = size(labels, ndims(labels));
        end

        % assign new parameter values now (don't wait for Solver.step)
        idx = [net.params.var];
        is_grad = ([net.params.trainMethod] == 1);
        net.setValue(idx(is_grad), w) ;

        % run network forward with new parameters
        net.eval({}, 'forward');

        % get value of last var (assumes only one output)
        pred = net.vars{end - 1};
        pred = reshape(pred, 1, 1, [], batch_size);  % reshape to 4D tensor for vl_nnloss

        % compute loss value
        switch o.loss
        case 'logistic'
          loss_value = vl_nnloss(pred, labels, 'loss', 'softmaxlog');
        otherwise
          error('Loss not yet supported with autolambda');
        end

        h_new = loss_value;

        % ratio between true curvature and quadratic fit curvature
        ratio = (h_new - h_curr) / M;

        % increase or decrease lambda based on ratio
        w1  = o.autolambda_w1 ^ o.autolambda_interval;
        thresh = o.autolambda_thresh;
        if isscalar(thresh)
          thresh = [thresh, 2 - thresh];
        end
        assert(diff(thresh) >= 0);

        if ratio < thresh(1)
          lambda = lambda / w1;
        elseif ratio > thresh(2)
          lambda = lambda * w1;
        end
        o.lambda = lambda;
        o.ratio = ratio;
      end

      o.count = o.count + 1;

    end

    function s = saveobj(o)
      % serialize to struct (called by the built-in function SAVE)
      % transfer state to CPU first
      s = o.saveGeneric();  % call parent class
      s.last_loss = gather(o.last_loss);
    end
  end

  methods (Static)
    function o = loadobj(s)
      % deserialize from struct (called by the built-in function LOAD)
      o = solvers.CurveBall();
      o.last_loss = s.last_loss;
      o.loadGeneric(s);  % call parent class

      % don't let parent class clear some parameters
      o.conserveMemory = false;
    end
  end
end

% computes the loss value, its gradient, and the hessian (multiplied by a
% vector x).
function [Hlx, Jl, loss_value] = hessian_grad_loss(net, loss, labels, x)
  if isvector(labels)  % categorical labels
    batch_size = numel(labels);
  else  % regression targets, get last dimension
    batch_size = size(labels, 4);
  end

  % get value of last var (assumes only one output)
  pred = net.vars{end - 1};
  pred = reshape(pred, 1, 1, [], batch_size);  % reshape to 4D tensor for vl_nnloss

  switch loss
  case 'ls'  % least-squares loss.
    % compute Hl * x. Hl = 2 / batch_size * I.
    Hlx = 2 / batch_size * x;

    % loss Jacobian, Jl
    assert(isequal(size(pred), size(labels)));
    Jl = 2 / batch_size * (pred - labels);

    % compute loss value
    loss_value = mean(sum((pred - labels).^2, 3), 4);

  case 'logistic'  % logistic loss.
    % compute Hl * x. for a single sample, Hl = diag(p) - p * p', where p
    % is a column-vector. for many samples, p has one column per sample,
    % and Hl is a block-diag with each block as above (so one independent
    % matrix-vec product per sample). first compute p' * x for all samples

    p = vl_nnsoftmax(pred);  % softmaxed probabilities

    if ismatrix(p)
      px = sum(p .* x, 1);
    else  % 4D tensor
      px = sum(p .* x, 3);
    end

    % now finish computing Hl * x = diag(p) * x - p * p' * x for every sample.
    Hlx = p .* x - p .* px;
    Hlx = Hlx / batch_size;

    if nargout <= 1, return, end  % early exit for single output

    % compute loss Jacobian, Jl. use backward pass of vl_nnloss
    Jl = vl_nnloss(pred, labels, single(1), 'loss', 'softmaxlog');  % loss gradient

    % compute loss value
    loss_value = vl_nnloss(pred, labels, 'loss', 'softmaxlog');

  otherwise
    error('Unknown loss.');
  end
end

repo: curveball


text:
function [rv,rvSS]=realized_threshold_multipower_variation(price,time,timeType,samplingType,samplingInterval,gamma,thresholdScale,thresholdType,subsamples)

c = thresholdScale;

logPrice =log(price);
% Filter prices and compute the RV
filteredLogPrice = realized_price_filter(logPrice,time,timeType,samplingType,samplingInterval);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filteredLogPrice  = log(price);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

returns = diff(filteredLogPrice);

L = 25;
V = inf*ones(size(returns));
K = -L:L;
K = 1/sqrt(2*pi)*exp(-(K/L).^2/2);
K(L+(-1:1)) = 0;
m = size(returns,1);
returns2 = returns.^2;

finished = false;
while ~finished
    ind = returns2<(c^2.*V);
    Vold  = V;
    for i=1:m
        pl = i-L:i+L;
        valid = pl>1 & pl<=m;
        tempInd = ind(pl(valid));
        w = K(valid);
        V(i) = w*(returns2(pl(valid)).*tempInd)/(w*tempInd);
    end
    if all(V==Vold)
        finished = true;
    end
end

expectedValue = 1./(2*normcdf(-c)*sqrt(pi))*(2/c^2) * gamma(3/2).*gammainc(c^2/2,3/2,'upper') * c^2 *V;

rv = returns(ind)'*returns(ind) + sum(expectedValue(~ind));

author: bashtage

repo: mfe-toolbox


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% A large number of Toolbox functions can operate on symbolic rather than
% numeric quantities.
%
% Consider the simple case of a rotation matrix

syms a
rotx(a)

% Or a more complex example for Euler angles

syms a b c
eul2r(a, b, c)

% Now let's consider a robot link with symbolic parameters

syms q A D alpha

L = Link('d', D, 'a', A, 'alpha', a, 'revolute');

% the link transform matrix for a joint angle q is then

L.A(q)

% Consider now a simple two-link robot, we load the model

mdl_twolink

% which has created a robot model in the workspace

twolink

% Now this is a numeric robot model, and we need to create a symbolic model

twolink_sym = twolink.sym()

% which appears very similar, however all the constants are now symbolic rather
% than numeric.

% Next define the two joint angles as symbolic variables

syms q1 q2

% and then the forward kinematics is

twolink_sym.fkine([q1, q2])

author: Allopart

repo: rbpf-gmapping


text:
function nlammap()
    %  NLAMMAP displays a map view of the seismicity in Lambert projection and ask for two input
    % points select with the cursor. These input points are
    % the endpoints of the crossection.
    %
    % Stefan Wiemer 2/95
    % updated: 12.10.2004, jochen.woessner@sed.ethz.ch
    % turned into function by Celso G Reyes 2017

    %global a
    %global main mainfault faults coastline vo s1_east s2_west s3_north s4_south
    %global mapl fipo
    %global h2 newa lat1 leng lon1 lon2 lat2

    %{
%% This is how you do it with the mapping toolbox -CGR
fig = figure;
axm=axesm('lambert','MapLatLimit',bounds2(ZG.primeCatalog.Latitude),...
  'MapLonLimit',bounds2(ZG.primeCatalog.Longitude));
plm=plotm(coastlat,coastlon,'k');
plm=plotm(ZG.primeCatalog.Latitude,ZG.primeCatalog.Longitude,'.');
disp('Select one end of cross-section')
p1=ginput(1); %TODO offer chance to redo/abort

disp('Select other end of cross-section')
p2=ginput(1);

    %}

    [c2, gcDist, zans] = plot_cross_section_from_mainmap; %was select_xsection();
    disp(c2)
    return

    ZG=ZmapGlobal.Data;
    report_this_filefun();
    %
    % Find out if figure already exists
    %
    mapl=findobj('Type','Figure','-and','Name','Seismicity Map (Lambert)');

    ZG.xsec_rotation_deg = 0;
    % Set up the Seismicity Map window Enviroment
    %
    if isempty(mapl)
        mapl = figure_w_normalized_uicontrolunits( ...
            'Name','Seismicity Map (Lambert)',...
            'NumberTitle','off', ...
            'backingstore','on',...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

        drawnow
        info2=@(s,ev)web(['file:' ZG.hodi '/zmapwww/chap4.htm#997433']);
        info1=@(s,ev)web(['file:' ZG.hodi '/zmapwww/chp11.htm#996756']) ;
        uicontrol('Units','normal',...
            'Position',[.0 .93 .08 .06],'String','Info ',...
            'callback',@info1_callback);

        uicontrol('Units','normal',...
            'Position',[.0 .93 .08 .06],'String','Info ',...
            'callback',@info2_callback);

    end % if figure exist
    %figure(mapl);
    %delete(findobj(mapl,'Type','axes'));
    plotmap(mapl);
    %{
    if isempty(coastline)
        coastline = [ZG.primeCatalog.Longitude(1) ZG.primeCatalog.Latitude(1)];
    end
    %}
    set(gca,'NextPlot','add')
    % Added try-catch to prevent failure if no coastline is inside
    % cross-section box, JW
    %try
    %{
    if length(coastline) > 1 %TODO what is coastline?
        lc_map(coastline(:,2),coastline(:,1),s3_north,s4_south,s1_east,s2_west)
        g = allchild(gca);
        set(g,'Color','k')

        %catch
    end
    set(gca,'NextPlot','add')
    try
        if length(faults) > 10
            lc_map(faults(:,2),faults(:,1),s3_north,s4_south,s1_east,s2_west)
        end
    catch
    end
    set(gca,'NextPlot','add')
    if ~isempty(mainfault)
        lc_map(mainfault(:,2),mainfault(:,1),s3_north,s4_south,s1_east,s2_west)
    end
    at_dep1 = ZG.primeCatalog.Depth<=dep1;
    at_dep2 = ZG.primeCatalog.Depth<=dep2 & ZG.primeCatalog.Depth>dep1;
    at_dep3 = ZG.primeCatalog.Depth<=dep3 & ZG.primeCatalog.Depth>dep2;
    if ZG.primeCatalog.Count > 5000
        %lc_event(ZG.primeCatalog.Latitude,ZG.primeCatalog.Longitude,'.k')
        lc_event(ZG.primeCatalog.Latitude(at_dep1),ZG.primeCatalog.Longitude(at_dep1),'.b',1);
        lc_event(ZG.primeCatalog.Latitude(at_dep2),ZG.primeCatalog.Longitude(at_dep2),'.g',1);
        lc_event(ZG.primeCatalog.Latitude(at_dep3),ZG.primeCatalog.Longitude(at_dep3),'.r',1);
    else
        lc_event(ZG.primeCatalog.Latitude(at_dep1),ZG.primeCatalog.Longitude(at_dep1),'+b');
        lc_event(ZG.primeCatalog.Latitude(at_dep2),ZG.primeCatalog.Longitude(at_dep2),'og');
        lc_event(ZG.primeCatalog.Latitude(at_dep3),ZG.primeCatalog.Longitude(at_dep3),'xr');

    end

    if ~isempty(ZG.maepi)
        lc_event(ZG.maepi.Latitude,ZG.maepi.Longitude,'hy',10,2.0)
    end
    if ~isempty(main)
        lc_event(main(:,2),main(:,1),'hk',10,2.0)
    end
    if ~isempty(vo)
        lc_event(vo.Latitude,vo.Longitude,'^r')
    end
    if ~isempty(well)
        lc_event(well(:,2),well(:,1),'dk')
    end
    %}
    labelList={'Select Endpoints by Mouse',...
        'Coordinate Input',...
        'Multiple segments',...
        'Rotate X-Section'};
    funList = { @do_nlammap2,... select endpoints by mouse
                @posinpu,...    coordinate input
                @musec,...      multiple segments
                @rotateit};     % rotate cross-section

    labelPos = [.05 .00 .40 .06];

    tmp1=ZG.primeCatalog.Latitude';
    tmp2=ZG.primeCatalog.Longitude';

    % dialog box for parameters
    zdlg=ZmapDialog();
    %zdlg.AddEdit(tag,label,value,tooltip);
    zdlg.AddEdit('xsec_width_km','Cross section width [km]',ZG.xsec_defaults.WidthKm,'cross section width, km');
    zdlg.AddPopup('uic',         'Selection Method:',labelList,1,...
        'Select an option for choosing cross section', funlist);
    zdlg.AddCheckbox('do_rotation','Rotate Cross Section', false, 'xsec_rotation_deg','Rotate Cross section');
    zdlg.AddEdit('xsec_rotation_deg','Rotation [deg]:',ZG.xsec_rotation_deg,'Rotate cross-section');
    [result,okPressed]=zdlg.Create('Name', 'Cross section parameters');
    if ~okPressed
        return
    end
    ZG.xsec_defaults.WidthKm = result.xsec_width_km;
    ZG.xsec_rotation_deg = result.xsec_rotation_deg;

    result.uic();

    function do_nlammap2()
        [xsecx xsecy,  inde] = mysect(tmp1,tmp2,ZG.primeCatalog.Depth,ZG.xsec_defaults.WidthKm);
        nlammap2;
    end

end

function rotateit()
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    if ~exist('lat1')
        errordlg('Create a cross-section first!');
    else
        [lat1, lat2, lon1, lon2] = rotate_xsection(lat1, lat2, lon1, lon2, ZG.xsec_rotation_deg);

        [xsecx xsecy,  inde] = mysect(tmp1,tmp2,ZG.primeCatalog.Depth,ZG.xsec_defaults.WidthKm,0,lat1,lon1,lat2,lon2);
        nlammap2;
    end

end

author: CelsoReyes


text:
function varargout = imarray(varargin)
% VL_IMARRAY  Flattens image array
%   J=VL_IMARRAY(A) creates a mosaic J from the array of images A. A
%   can be either a M*N*K array, storing one gray-scale image per
%   slice, or a M*N*3*K or M*N*K*3 array, storing a true color RGB
%   image per slice. The function returns an image J which is a tiling
%   of the images in the array. Tiles are filled from left to right
%   and top to bottom.
%
%   VL_IMARRAY(...) displays the image J rather than returning it.
%
%   VL_IMARRAY() accepts the following options:
%
%   Spacing:: 0
%     Separate the images by a border of the specified width (the
%     border is assigned 0 value, which usually corresponds to black).
%
%   Layout:: empty
%     Specify a vector [TM TN] with the number of rows and columns of
%     the tiling. If equal to [] the layout is computed automatically.
%
%   Movie:: false
%     Display or return a movie instead of generating a tiling.
%
%   CMap:: []
%     Specify a colormap to construct a movie when the input is an
%     indexed image array. If not specified, MATLAB default colormap
%     is used.
%
%   Reverse:: true
%     Start filling the mosaic tiles from the bottom rather than from
%     the top.
%
%   See also:VL_IMARRAYSC(), VL_HELP().
[varargout{1:nargout}] = vl_imarray(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function halton_test09 ( )

%*****************************************************************************80
%
%% TEST09 tests SPHERE_UNIT_HALTON_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num2 = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST09\n' );
  fprintf ( 1, '  For the unit sphere in 2 dimensions (the circle):\n' );
  fprintf ( 1, '  HALTON generates "U1" points,\n' );
  fprintf ( 1, '  U1_TO_SPHERE_UNIT_2D samples the circle;\n' );

  dim_num = 1;
  halton_dim_num_set ( dim_num );
  n = 5;
  step = 0;
  halton_step_set ( step );
  seed(1:dim_num) = 0;
  halton_seed_set ( seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  DIM_NUM = %12d\n', dim_num );
  fprintf ( 1, '  N =    %12d\n', n );
  fprintf ( 1, '  STEP = %12d\n', step );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A few sample values:\n' );
  fprintf ( 1, '\n' );

  for j = 1 : n
    u = halton ( );
    x = u1_to_sphere_unit_2d ( u );
    fprintf ( 1, '  ' );
    for i = 1 : dim_num2
      fprintf ( 1, '%6f  ', x(i) );
    end
    fprintf ( 1, '\n' );
  end

  n = 1000;
  step = 0;
  halton_step_set ( step );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N =    %12d\n', n );
  fprintf ( 1, '  STEP = %12d\n', step );

  average(1:dim_num2) = 0.0;

  for i = 1 : n
    u = halton ( );
    x = u1_to_sphere_unit_2d ( u );
    average(1:dim_num) = average(1:dim_num) + x(1:dim_num);
  end

  average(1:dim_num2) = average(1:dim_num2) / n;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average the points, which should get a value\n' );
  fprintf ( 1, '  close to zero, and closer as N increases.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average:        ' );
  for i = 1 : dim_num2
    fprintf ( 1, '  %12f', average(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now choose a random direction, sample the same\n' );
  fprintf ( 1, '  number of points, and compute the dot product with\n' );
  fprintf ( 1, '  the direction.\n' );
  fprintf ( 1, '  Take the absolute value of each dot product\n' );
  fprintf ( 1, '  and sum and average.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We expect a value near 2 / PI = 0.6366...\n' );

  for j2 = 1 : 5

    seed_val = get_seed ( );
    step = seed_val + 111 * j2;
    halton_step_set ( step );

    u = halton ( );
    v = u1_to_sphere_unit_2d ( u );

    step = 0;
    halton_step_set ( step );

    dot_average = 0.0;

    for j = 1 : n
      u = halton ( );
      x = u1_to_sphere_unit_2d ( u );
      dot_average = dot_average + abs ( x(1:dim_num2) * v(1:dim_num2)' );
    end

    dot_average = dot_average / n;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Random V:         ' );
    for i = 1 : dim_num2
      fprintf ( 1, '  %12f', v(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Average |(XdotV)| = %f\n', dot_average );

  end

  return
end