text:
%   CLASS Galileo_SS
% =========================================================================
%
% DESCRIPTION
%   container of Galileo Satellite System parameters
%
% REFERENCES
%   CRS parameters, according to each GNSS system CRS definition
%   (ICD document in brackets):
%
%   *_GAL --> GTRF     (Galileo-ICD 1.1)
%   Standard: https://www.gsc-europa.eu/system/files/galileo_documents/Galileo_OS_SIS_ICD.pdf
%
%   Other useful links
%     - http://www.navipedia.net/index.php/Galileo_Signal_Plan
%     - http://www.navipedia.net/index.php/Reference_Frames_in_GNSS
%    Ellipsoid definition is actually coming from this presentation:
%     - http://gage6.upc.es/eknot/Professional_Training/PDF/Reference_Systems.pdf
%    at the moment of writing the Galileo Geodetic Reference Service Provider (GRSP) website (http://www.ggsp.eu/)
%    is actually offline, and the GTRF16v01 (or any other RF) cannot be found online.
%    Since GTRF seems to be based on ITRF (that does not define a reference ellispoid) http://itrf.ign.fr/faq.php?type=answer
%    The ellipsoid found in the presentation is used as reference
%
%   Useful data for the future (PPP):
%     - https://www.gsc-europa.eu/support-to-developers/galileo-iov-satellite-metadata#3.2

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

classdef Galileo_SS < Satellite_System
    properties (Constant, Access = 'public')
        SYS_EXT_NAME = 'Galileo'; % full name of the constellation
        SYS_NAME     = 'GAL';     % 3 characters name of the constellation, this "short name" is used as fields of the property list (struct) to identify a constellation
        SYS_C        = 'E';       % Satellite system (ss) character id

        % System frequencies as struct [MHz]
        F = struct('E1', 1575.420, ...
                   'E5a', 1176.450, ...
                   'E5b', 1207.140, ...
                   'E5', 1191.795, ...
                   'E6', 1278.750)

        % Array of supported frequencies [MHz]
        F_VEC = struct2array(Galileo_SS.F) * 1e6;

        % Array of the corresponding wavelength - lambda => wavelengths
        L_VEC = 299792458 ./ Galileo_SS.F_VEC;

        N_SAT = 36;       % Maximum number of satellite in the constellation
        PRN = (1 : 36)';  % Satellites id numbers as defined in the constellation

        % CODE2DATA ftp://igs.org/pub/data/format/rinex303.pdf
        CODE_RIN3_ATTRIB  = {'ZXACB F' 'XQI F' 'XQI F', 'XQI F', 'ZXACB F'}; % last letter of the observation code Assumption: Public regualted service (PRS ) better than Pilot (C) better than data (A), pilot channel seems to be the quadra pahse one
        CODE_RIN3_DEFAULT_ATTRIB  = {'C' 'I' 'I' 'I' 'C'}; % last letter of the observation code
        CODE_RIN3_2BAND  = '15786';             % id for the freq as stored in F_VEC
        IONO_FREE_PREF  = ['15'; '18'; '17'; '16'; '56'; '76'; '86'; '57'; '58'; '78'];  % to be evaluated which combination is really better
    end

    properties (Constant, Access = 'private')
        % GPS (WGS84) Ellipsoid semi-major axis [m]
        ELL_A = 6378137;
        % GPS (WGS84) Ellipsoid flattening
        ELL_F = 1/298.257222101;
        % GPS (WGS84) Ellipsoid Eccentricity^2
        ELL_E2 = (1 - (1 - Galileo_SS.ELL_F) ^ 2);
        % GPS (WGS84) Ellipsoid Eccentricity
        ELL_E = sqrt(Galileo_SS.ELL_E2);
    end

    properties (Constant, Access = 'public')
        % Structure of orbital parameters (ellipsoid, GM, OMEGA_EARTH_DOT)
        ORBITAL_P = struct('GM', 3.986004418e14, ...               % Galileo (Galileo-ICD) Gravitational constant * (mass of Earth) [m^3/s^2]
                                    'OMEGAE_DOT', 7.2921151467e-5, ...      % Galileo (Galileo-ICD) Angular velocity of the Earth rotation [rad/s]
                                    'ELL',struct( ...                       % Ellipsoidal parameters Galileo (GTRF)
                                    'A', Galileo_SS.ELL_A, ...          % Ellipsoid semi-major axis [m]
                                    'F', Galileo_SS.ELL_F, ...          % Ellipsoid flattening
                                    'E', Galileo_SS.ELL_E, ...          % Eccentricity
                                    'E2', Galileo_SS.ELL_E2));          % Eccentricity^2
        ORBITAL_INC = 56;    % Orbital inclination
        ORBITAL_RADIUS  = 23222000 + 6378137; % Orbital radius
    end

    methods
        function this = Galileo_SS(offset)
            % Creator
            % SYNTAX: Galileo_SS(<offset>);
            if (nargin == 0)
                offset = 0;
            end
            this@Satellite_System(offset);
        end

        function copy = getCopy(this)
            % Get a copy of this
            copy = Galileo_SS(this.getOffset());
            copy.import(this);
        end

    end
end

author: goGPS-Project


text:
function linplus_test195 ( )

%*****************************************************************************80
%
%% TEST195 tests R8CI_EVAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST195\n' );
  fprintf ( 1, '  R8CI_EVAL finds the eigenvalues of\n' );
  fprintf ( 1, '  a real circulant system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r8ci_random ( n, seed );

  r8ci_print ( n, a, '  The circulant matrix:' );

  lambda = r8ci_eval ( n, a );

  c8vec_print ( n, lambda, '  The eigenvalues:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ a ] = e_qgreedy(Q,s,epsilon )
%E_QGREEDY
%
% input -------------------------------------------------------------------
%
%       o Q       : (num_states x num_actions), Q-value table
%
%       o s       : (1 x 1), current state index
%
%       o epsilon : (1 x 1), exploration noise [0,1]
%

actions = size(Q,2);
if (rand()>epsilon)
    [~,a] = max(Q(s,:));
else
    a = randi(actions);
end

end

author: epfl-lasa


text:
function varargout = run_opf(varargin)
%run_opf  Run an optimal power flow.
%
%   run_opf(d, mpopt)
%   run_opf(d, mpopt, ...)
%   task = run_opf(...)
%
%   See also run_mp.

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

[varargout{1:nargout}] = run_mp(@mp.task_opf, varargin{:});

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%CONE Collision checking class for cone-type primitive shapes
%
% Subclass of the Curvilinear class. The base of the cone intersects
% the transformation frame's origin. May be circular or elliptical.
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% Syntax:
%  (1) cone = Cylinder(T, s)
%  (2) cone = Cone(T, s, ...)
%  (3) cone = Cone(P, r, ...)
%  (3) cone = Cone()
%
%  (2) is as per (1) but with other properties set in name-value pairs
%  (3) is as per (1) or (2) but is defined by two points and a radius -
%       the cross-section thus must be circular and not elliptical
%  (4) returns a cone with default values
%
% Outputs:
%  cone : Cone object
%
% Inputs:
%  T      : Transformation matrix. If empty, defaults to eye(4).
%  s      : Scale vector. If empty, defaults to [1 1 1].
%  P      : Cell, whose first element is the origin, and the second
%            element is the other endpoint. Cell and points may be
%            either column or row vectors.
%  r      : Scalar radius value, for syntax mode (3)
%  ...    : Options - other properties in name-value pairs
%
% See documentation for information on properties and methods
% (type doc Shape into the command window)
%
% See also Box.Box Curvilinear Cylinder.Cylinder Ellipsoid.Ellipsoid
% Shape Sphere.Sphere

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

classdef Cone < Curvilinear

    methods
        function s = Cone(T, S, varargin)
            if ~nargin
               T = [];
               S = [];
            end

            p = @(t)(1-t);

            coneDefaults = [{'n',2,'faces',[1 0]}, varargin];
            s = s@Curvilinear(T,S,p,coneDefaults{:});
        end
    end
end

author: Allopart

repo: rbpf-gmapping


text:
function UpdateClustersViewGUI(hfig)
clusID = getappdata(hfig,'clusID');
clusgroupID_view = getappdata(hfig,'clusgroupID_view');
i_fish = getappdata(hfig,'i_fish');
global VAR;
ClusGroup = VAR(i_fish).ClusGroup{clusgroupID_view};

global hclusname hclusmenu;
set(hclusname,'String',ClusGroup(clusID).name);
menu = MakeNumberedMenu({ClusGroup.name});
set(hclusmenu,'String', menu,'Value',clusID+1);
end

author: xiuyechen

repo: FishExplorer


text:
function p = spm_LAP_pg(x,v,h,M);
% default precision function for LAP models (hidden states)
% FORMAT p = spm_LAP_pg(x,v,h,M);
%
% x  - hidden states
% v  - causal states
% h  - precision parameters
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_LAP_pg.m 3694 2010-01-22 14:16:51Z karl $

% fixed components
%--------------------------------------------------------------------------
p = sparse(M.n,1);
try
    W = diag(M.W);
    if all(W)
        p = log(W);
    end
end

% free components
%--------------------------------------------------------------------------
for i = 1:length(M.R)
    p = p + h(i)*diag(M.R{i});
end

author: spm


text:
function [ matlabdatenum ] = rawjavacalendar2datenum( cal )
%   Converts a java.util.Calendar date local time into a Matlab datenum.
%   Keeps the original local time, does not try to convert to UTC.
    javaSerialDate = cal.getTimeInMillis() + cal.get(cal.ZONE_OFFSET) + cal.get(cal.DST_OFFSET);
    matlabdatenum = datenum([1970 1 1 0 0 javaSerialDate / 1000]);
end

%http://www.mathworks.co.uk/support/solutions/en/data/1-9B9H2S/index.html?product=ML&solution=1-9B9H2S

author: Sable

repo: mcbench-benchmarks


text:
% rf_addGrad.m
% Jamie Near, McGill University 2020.
%
% USAGE:
% RF_out=rf_addGrad(RF_in,grad);
%
% DESCRIPTION:
% Add a gradient waveform to the input RF pulse.
%
% INPUTS:
% RF_in     = Input RF pulse definition structure
% grad      = If grad is a scalar, then rf_addGrad will add a constant
%             gradient with amplitude equal to the value of grad.  If grad
%             is a vector, then rf_addGrad will add the specified gradient
%             vector to the RF waveform.  In either case, the gradient
%             should be specified in units of [G/cm].
%
% OUTPUTS:
% RF_out    = Output rf waveform following the addition of gradient
%             waveform.

function RF_out=rf_addGrad(RF_in,grad);

if ~isstruct(RF_in)
    error('ERROR:  the input RF pulse must be in structure format.  Try using rf_readwaveform to convert it!  Aborting.  ');
end

newWaveform=RF_in.waveform;

%Check if the input RF pulse already has a gradient.
if size(newWaveform,2)>3
    keepGoing='n';
    disp('WARNING:  Input waveform already has a gradient waveform!');
    keepGoing=input('Do you wish to overwrite the existing gradient waveform (y or [n])','s');

    if strcmp(keepGoing,'y') || strcmp(keepGoing,'Y');
        disp('OK.  Overwriting existing gradient waveform.');
    elseif strcmp(keepGoing,'n') || strcmp(keepGoing,'N');
        error('ABORTING!!!');
    else
        error('Response not recognized.  ABORTING!!');
    end
end

%Now add the gradient.
if isscalar(grad)
    %Adding a scalar gradient
    newWaveform(:,4)=grad*ones(size(newWaveform,1),1);
else
    %Adding a gradient waveform
    if size(grad) ~= size(newWaveform,1)
        if size(grad') ~= size(newWaveform,1);
            error('ERROR:  Gradient waveform does not match the length of the input RF pulse waveform!!  ABORTING!!');
        end
        newWaveform(:,4)=grad';
    end
    newWaveform(:,4)=grad';
end

%now it's time to find out the time-bandwidth product:
%First make a high resolution plot the pulse profile over a wide bandwidth:
[mv,sc]=bes(rf,Tp*1000,'f',w1max/1000,-5+f0/1000,5+f0/1000,100000);
if isstr(RF_in.type)
    if RF_in.type=='exc'
        index=find(mv(3,:)<0.5);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    elseif RF_in.type=='ref'
        index=find(mv(3,:)<0);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    elseif RF_in.type=='inv'
        index=find(mv(3,:)<0);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    end
elseif isnumeric(RF_in.type)
    mz=cos(type); %Find out the Mz value immediately following the pulse.
    thr=(1+mz)/2;  %Find out the Mz value mid-way between 1 and the mz (half-max):
    index=find(mv(3,:)<thr);  %Find the indices of the corresponding "full width"
    bw=sc(index(end))-sc(index(1));  %Now find the bandwidth at that point ("Full width at half max").
end

%Now make a very high resolution plot the pulse profile over a narrower bandwidth:
[mv,sc]=bes(rf,Tp*1000,'f',w1max/1000,-bw+f0/1000,bw+f0/1000,100000);
if isstr(RF_in.type)
    if RF_in.type=='exc'
        index=find(mv(3,:)<0.5);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    elseif RF_in.type=='ref'
        index=find(mv(3,:)<0);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    elseif RF_in.type=='inv'
        index=find(mv(3,:)<0);
        bw=sc(index(end))-sc(index(1));
        %plot(sc(index),mv(3,index),'.-',sc,mv(3,:));
    end
elseif isnumeric(RF_in.type)
    mz=cos(RF_in.type); %Find out the Mz value immediately following the pulse.
    thr=(1+mz)/2;  %Find out the Mz value mid-way between 1 and the mz (half-max):
    index=find(mv(3,:)<thr);  %Find the indices of the corresponding "full width"
    bw=sc(index(end))-sc(index(1));  %Now find the bandwidth at that point ("Full width at half max").
end

RF_out=RF_in;
RF_out.waveform=newWaveform;
RF_out.isGM=true;
RF_out.tbw='N/A - gradient modulated pulse';
RF_out.isGM=true;
RF_out.tthk=bw*Tp; %This is the time x sliceThickness product for
                   %gradient modulated pulses.  It is in units [cm.s]

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function read_ceos_raw ( output_file_pre, start_line_blk, blk )
%
% -------------------------------------------------------------------------
% This function reads the CEOS format RSAT-1 RAW data from the CD
% This function is called by extract_data.m
% -------------------------------------------------------------------------
% output_file_pre : output file prefix ( path + prefix + block number)
% start_line_blk  : extract data starting from row number 'start_line_blk'
% blk             : the block number
%
% input_path     : Input file path
% length_replica : Length (I & Q) of replica (for Fine Beam -> 2880 bytes)
% tot_Nrg_cells  : Total number of columns available in the dataset
% Nrg_lines_blk  : Number of samples in azimuth
% first_rg_cell  : Extract data starting from column number 'first_rg_cell'
% Nrg_cells      : Number of samples in ranged
% UseMATfiles    : Do you want to write MAT files for the data rather than
%                  binary files for all the data and parameters?
% -------------------------------------------------------------------------

% date created  : May 05, 2004  by  Kaan Ersahin & Millie Sikdar
% date modified : Nov 30, 2004  by  Kaan Ersahin
% Modified:  Dec 17 by Ian Cumming
%             - Added UseMATfiles option
% -------------------------------------------------------------------------

% 'RSI - D4' refers to the document RSI-GS-026, Rev 3.0 - May 8, 2000
%  that describes the detailed format of the RADARSAT-1 raw data CDs.
tic
num_aux            = 50;                        % RSI - D4 Pg.32
num_header         = 192;                       % RSI - D4 Pg.32
file_header_length = 16252;                     % RSI - D4 Pg.103
load CD_run_params

%  Specify input data file path and open to read
file_in  = strcat( input_path, CD_data_file_name );
fid      = fopen( file_in, 'r' );

%  Create output file names and open to write if MAT files not used
if ~UseMATfiles
    file_replica 	 = '_replica.dat';     % one replica per 8-line
    file_data_header = '_data_header.dat'; % header of each line
    file_aux     	 = '_aux.dat';         % auxiliary data from each line
    file_data    	 = '_data.dat';        % radar signal data
    file_file_header = '_file_header.dat'; % header of file 'dat_01.001'

    file_replica     = strcat(output_file_pre, file_replica);
    file_data_header = strcat(output_file_pre, file_data_header);
    file_aux         = strcat(output_file_pre, file_aux);
    file_data        = strcat(output_file_pre, file_data);
    file_file_header = strcat(output_file_pre, file_file_header);

    fod1             = fopen( file_replica    , 'w' );
    fod2             = fopen( file_data_header, 'w' );
    fod3             = fopen( file_aux        , 'w' );
    fod4             = fopen( file_data       , 'w' );
    fod5             = fopen( file_file_header, 'w' );

    %  Read file header
    file_header      = fread( fid, [file_header_length,1], 'uint8' );
    %  Write file header in a binary file as unsigned integers
    fwrite( fod5, file_header, 'uint8');
    fclose( fod5 );
end  % of if ~UseMATfiles

% some useful calculations
num_pixel_data   = 2*tot_Nrg_cells;                               %  18576
h                = num_header;                                    %    192
ha               = num_header + num_aux;                          %    242
hpa              = num_header + num_aux + num_pixel_data;         %  18818

% Define 'rep_block' as the 8-line block of data including one replica
% 'start_rep' and 'end_rep' locates the replica in the 'rep_block'
% 'start_col' and 'end_col' locates the area of interest in the columns
rep_block_length = length_replica + 8*hpa ;

start_rep        = (first_replica-1)*hpa + num_header + num_aux + 1;
end_rep          = start_rep + length_replica - 1 ;

start_col        = ha + 2*(first_rg_cell - 1) + 1 ;
end_col          = start_col + 2*Nrg_cells - 1;

% The number of bytes to skip from the beginning of the file (Nbytes_skip)
Nrep_blks_row1   = ceil(start_line_blk/8) - 1;
Nbytes_blocks    = Nrep_blks_row1 * rep_block_length;
Nbytes_skip      = file_header_length + Nbytes_blocks;

fseek( fid, Nbytes_skip , 'bof' );  % go to start of the 1st data block
data = int8( zeros(Nrg_lines_blk, 2*Nrg_cells) );   % Allocate data array

N_8line_blocks = ceil(Nrg_lines_blk/8);
fprintf('\nReading %1.0f small 8-line blocks from range line %5.0f\n',...
    N_8line_blocks, start_line_blk )

for kb = 1 : N_8line_blocks   %  Read and decode 8 lines at a time

    % read one 8-line block of radar data and replica
    temp = uint8( fread( fid, [rep_block_length,1], 'uint8' ) );

    % separate the replica and 8 lines of [header, aux and data]
    replica = temp(start_rep : end_rep)';
    temp = [ temp(1 : start_rep-1)' temp(end_rep+1 : rep_block_length)' ];
    temp = reshape(temp, hpa, 8)';   % Shape array as 8 rows by all columns

    % Extract the desired set of range cells
    data(1+8*(kb-1):8*kb, :) = temp( : , start_col : end_col );

    if ~UseMATfiles  % Write replica, header, aux and data in a binary file
        header (:,1:8) = temp( : ,         1 : h       )';
        aux    (:,1:8) = temp( : ,       h+1 : ha      )';
        count1 = fwrite( fod1, replica, 'ubit4' ); % 4-bit replica
        count2 = fwrite( fod2, header , 'uint8' ); % 8-bit header
        count3 = fwrite( fod3, aux    , 'uint8' ); % 8-bit aux data
        count4 = fwrite( fod4, data(1+8*(kb-1):8*kb, :)'   , 'ubit4' );
    end  % of if UseMATfiles
end

fclose(fid);           % Close input file
if ~UseMATfiles        % Close output files
    fclose(fod1);   fclose(fod2);   fclose(fod3);   fclose(fod4);
    fprintf('This block written to binary file:\n%s\n', file_data )
end

%--------------------------------------------------------------------------
%  Decode the data, save as MAT file , if UseMATfiles == 1
%--------------------------------------------------------------------------

if UseMATfiles == 1
    disp 'Decode data by subtracting 16 from upper register'
    vv    = version;   vers = str2num(vv(1));  % Find MATLAB version
    seven = int8( 7 );

    if vers == 7
        one     = int8( 1 );
        two     = int8( 2 );
        sixteen = int8( 16 );
        data    = two*(data - sixteen*int8(data > seven)) + one;
    else
        data    = int8( 2*(double(data) - 16*double(data > seven) ) + 1 );
    end

    %  Separate the I and Q channels and form a complex array
    data = complex( data(:, 1:2:2*Nrg_cells), data(:, 2:2:2*Nrg_cells) );

    if SaveV6
       eval(['save -v6 ' output_path 'CDdata' num2str(blk) ' data'])
       disp  'Save int8 data in MAT file in V6 format'
    else
       eval(['save ' output_path 'CDdata' num2str(blk) ' data'])
       fprintf('Save int8 data in file:  CDdata%1.0f.mat\n', blk)
    end
    % somedata = data(1,1:4)  % Check a few data samples
end   % of if UseMATfiles
toc

author: denkywu

repo: SAR-Synthetic-Aperture-Radar


text:
function [net, info] = cnn_mnist(varargin)
% CNN_MNIST  Demonstrated MatConNet on MNIST

run(fullfile(fileparts(mfilename('fullpath')),...
  '..', 'matlab', 'vl_setupnn.m')) ;

opts.expDir = fullfile('data','mnist-baseline') ;
opts.gpus = [];
opts.pushbullet = []; % usage : pbNotify('accessToken'); from https://www.pushbullet.com/account
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.dataDir = fullfile('data','mnist') ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.useBnorm = false ;
opts.train.batchSize = 100 ;
opts.train.numEpochs = 20 ;
opts.train.continue = true ;
opts.train.gpus = opts.gpus;
opts.train.learningRate = 0.001 ;
opts.train.expDir = opts.expDir ;
opts.train.pushbullet = opts.pushbullet ;
opts = vl_argparse(opts, varargin) ;

% --------------------------------------------------------------------
%                                                         Prepare data
% --------------------------------------------------------------------

if exist(opts.imdbPath, 'file')
  imdb = load(opts.imdbPath) ;
else
  imdb = getMnistImdb(opts) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

net = cnn_mnist_init('useBnorm', opts.useBnorm) ;

% --------------------------------------------------------------------
%                                                                Train
% --------------------------------------------------------------------
% b = [10 10] + [10 10 10] ; error test
[net, info] = cnn_train(net, imdb, @getBatch, ...
    opts.train, ...
    'val', find(imdb.images.set == 3)) ;

if ~isempty(opts.pushbullet)
    opts.pushbullet.notify('training complete') ;
end
% --------------------------------------------------------------------
function [im, labels] = getBatch(imdb, batch)
% --------------------------------------------------------------------
im = imdb.images.data(:,:,:,batch) ;
labels = imdb.images.labels(1,batch) ;

% --------------------------------------------------------------------
function imdb = getMnistImdb(opts)
% --------------------------------------------------------------------
% Preapre the imdb structure, returns image data with mean image subtracted
files = {'train-images-idx3-ubyte', ...
         'train-labels-idx1-ubyte', ...
         't10k-images-idx3-ubyte', ...
         't10k-labels-idx1-ubyte'} ;

if ~exist(opts.dataDir, 'dir')
  mkdir(opts.dataDir) ;
end

for i=1:4
  if ~exist(fullfile(opts.dataDir, files{i}), 'file')
    url = sprintf('http://yann.lecun.com/exdb/mnist/%s.gz',files{i}) ;
    fprintf('downloading %s\n', url) ;
    gunzip(url, opts.dataDir) ;
  end
end

f=fopen(fullfile(opts.dataDir, 'train-images-idx3-ubyte'),'r') ;
x1=fread(f,inf,'uint8');
fclose(f) ;
x1=permute(reshape(x1(17:end),28,28,60e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 't10k-images-idx3-ubyte'),'r') ;
x2=fread(f,inf,'uint8');
fclose(f) ;
x2=permute(reshape(x2(17:end),28,28,10e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 'train-labels-idx1-ubyte'),'r') ;
y1=fread(f,inf,'uint8');
fclose(f) ;
y1=double(y1(9:end)')+1 ;

f=fopen(fullfile(opts.dataDir, 't10k-labels-idx1-ubyte'),'r') ;
y2=fread(f,inf,'uint8');
fclose(f) ;
y2=double(y2(9:end)')+1 ;

set = [ones(1,numel(y1)) 3*ones(1,numel(y2))];
data = single(reshape(cat(3, x1, x2),28,28,1,[]));
dataMean = mean(data(:,:,:,set == 1), 4);
data = bsxfun(@minus, data, dataMean) ;

imdb.images.data = data ;
imdb.images.data_mean = dataMean;
imdb.images.labels = cat(2, y1, y2) ;
imdb.images.set = set ;
imdb.meta.sets = {'train', 'val', 'test'} ;
imdb.meta.classes = arrayfun(@(x)sprintf('%d',x),0:9,'uniformoutput',false) ;

author: thomas-koehler

repo: SupER


text:
% Mathematics Q2410741
% https://math.stackexchange.com/questions/2410741
% Matrix differentiation on ∥A−B∘X∥2F+λ∥X∥2F
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     30/08/2017
%   *   First release.

%% General Parameters

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = OFF;

DIFF_MODE_FORWARD   = 1;
DIFF_MODE_BACKWARD  = 2;
DIFF_MODE_CENTRAL   = 3;

%% Simulation Parameters

numRows = 2;

paramLambda = 0.1 * rand([1, 1]);

difMode = DIFF_MODE_FORWARD;
epsVal  = 1e-6;

numSamples      = 1000;
numIterations   = 2500;
stepSize        = 0.001;

%% Generate Data

vX = rand([numRows, 1]);

mA = randn([numRows, numRows]);
mB = randn([numRows, numRows]);

hObjFun = @(vX) (norm(mA - (mB .* (vX * vX.')), 'fro') .^ 2) + (paramLambda * (norm(vX * vX.', 'fro') ^ 2));
% hObjFun = @(vX) sum(sum((mA - mB .* (vX * vX.')) .^ 2)) + (paramLambda * sum(sum((vX * vX.') .^ 2)));

%% Numerical Solution

vGNumerical = CalcFunGrad(vX, hObjFun, difMode, epsVal);

%% Analytic Solution

vG = zeros([numRows, 1]);

for kk = 1:numRows
    for ii = 1:numRows
        for jj = 1:numRows
            if(ii == kk)
                vG(kk) = vG(kk) + (-2 * mB(ii, jj) * vX(jj) * (mA(ii, jj) - mB(ii, jj) * vX(ii) * vX(jj)));
                vG(kk) = vG(kk) + (2 * paramLambda * vX(jj) * (vX(ii) * vX(jj)));
            end
            if(jj == kk)
                vG(kk) = vG(kk) + (-2 * mB(ii, jj) * vX(ii) * (mA(ii, jj) - mB(ii, jj) * vX(ii) * vX(jj)));
                vG(kk) = vG(kk) + (2 * paramLambda * vX(ii) * (vX(ii) * vX(jj)));
            end
        end
    end
end

%% Analysis

disp(['Maximum Error Between Numercial Differntiation to Analytic - ', num2str(max(abs(vG - vGNumerical)))]);

%% Minimization

mZ      = zeros([numSamples, numSamples]);
vXGrid  = linspace(-5, 5, numSamples);

for jj = 1:numSamples
    for ii = 1:numSamples
        mZ(ii, jj) = hObjFun([vXGrid(jj); vXGrid(ii)]);
    end
end

[minVal, argMin] = min(mZ(:));
[minValRowIdx, minValColIdx] = ind2sub([numSamples, numSamples], argMin);

% Numercial Optimization

mX          = zeros([numRows, numIterations]);
vObjFunVal  = zeros([numIterations, 1]);

% Initialization
vX              = rand([numRows, 1]);
vX              = [-2; 4];

mX(:, 1)        = vX;
vObjFunVal(1)   = hObjFun(vX);

for ii = 2:numIterations
    vG              = CalcFunGrad(vX, hObjFun, difMode, epsVal);
    vX              = vX - (stepSize * vG);
    mX(:, ii)       = vX;
    vObjFunVal(ii)  = hObjFun(vX);
end

figureIdx = figureIdx + 1;

hFigure         = figure('Position', [100, 100, 720, 960]);
hAxes           = subplot(3, 1, 1:2);
set(hAxes, 'NextPlot', 'add');
hImageObject    = image('XData', vXGrid, 'YData', vXGrid, 'CData', mZ, 'CDataMapping', 'scaled');
[mCountourLines, hContouObj]  = contour(vXGrid, vXGrid, mZ);
set(hContouObj, 'LineColor', 'white');
hLineSeries(1)  = line(vXGrid(minValColIdx), vXGrid(minValRowIdx));
set(hLineSeries(1), 'LineStyle', 'none', 'Marker', 'diamond', 'MarkerSize', markerSizeLarge, 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'r');
hLineSeries(2)  = line(mX(2, :), mX(1, :));
set(hLineSeries(2), 'LineStyle', 'none', 'Marker', 'o', 'MarkerSize', 5, 'MarkerEdgeColor', 'm', 'MarkerFaceColor', 'w');
set(hAxes, 'DataAspectRatio', [1, 1, 1]);
set(get(hAxes, 'Title'), 'String', {['Objective Function Value']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['x_2']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['x_1']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend([hLineSeries(1), hLineSeries(2)], {['Optimal Argument'], ['Gradient Descent Iterations']});

hAxes = subplot(3, 1, 3);
hLineSeries = plot(1:numIterations, vObjFunVal);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['Objective Function Value']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Iteration Number']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Value']}, ...
    'FontSize', fontSizeAxis);

saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
classdef EdgesConnectivitiesComputer < handle

    properties (GetAccess = public, SetAccess = private)
        nodesInEdges
        edgesInElem
        nEdgeByElem
        nNodeByEdge
        localNodeByEdgeByElem
    end

    properties (Access = private)
        allToUnique
        uniqueToAll
        nodesInAllEdges
    end

    properties (Access = private)
        nodesByElem
        nElem
        nAllEdges
        localEdgesInElem
    end

    methods (Access = public)

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

        function compute(obj)
            obj.computeNodesInAllEdges();
            obj.computeUniqueEdgesIndex();
            obj.computeNodesInUniqueEdges();
            obj.computeEdgesInElem();
            obj.computeLocalOrientedEdgeConnec();
        end

        function cV = computeConnectedVertex(obj,vertex)
            e  = obj.computeAllEdgesOfVertex(vertex);
            cV = obj.computeOtherVertexOfEdge(e,vertex);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.nodesByElem = cParams.nodesByElem;
            nNodes = size(cParams.nodesByElem,2);
            switch nNodes
                case 3
                    obj.localEdgesInElem = nchoosek(1:nNodes,2);%[1 2; 2 3; 3 1];
                    obj.localEdgesInElem = [1 2; 2 3; 3 1];
                case 4
                    obj.localEdgesInElem =  nchoosek(1:nNodes,2);%[1 2; 2 3; 3 1];
            end
            obj.nElem = size(obj.nodesByElem,1);
            obj.nEdgeByElem = size(obj.localEdgesInElem,1);
            obj.nNodeByEdge = size(obj.localEdgesInElem,2);
            obj.nAllEdges = obj.nElem*obj.nEdgeByElem;
        end

        function computeNodesInAllEdges(obj)
            nodesInE = zeros(obj.nAllEdges,obj.nNodeByEdge);
            allElem = 1:obj.nElem;
            for inode = 1:obj.nNodeByEdge
                for iedge = 1:obj.nEdgeByElem
                    lNode = obj.localEdgesInElem(iedge,inode);
                    index = (allElem - 1)*obj.nEdgeByElem + iedge;
                    nodes = obj.nodesByElem(:,lNode);
                    nodesInE(index,inode) = nodes;
                end
            end
            obj.nodesInAllEdges = nodesInE;
        end

        function computeUniqueEdgesIndex(obj)
            nE = obj.nodesInAllEdges;
            [nE] = sort(nE,2);
            [~,a2u,u2a] = unique(nE,'rows');
            obj.allToUnique = a2u;
            obj.uniqueToAll = u2a;
        end

        function computeNodesInUniqueEdges(obj)
            a2u = obj.allToUnique;
            nodesInUniqueEdges = obj.nodesInAllEdges(a2u,:);
            obj.nodesInEdges = nodesInUniqueEdges;
        end

        function computeEdgesInElem(obj)
            uniqueEdgeByElem = obj.computeUniqueEdgesByElem();
            eInElem = zeros(obj.nElem,obj.nEdgeByElem);
            for iEdge = 1:obj.nEdgeByElem
                edge = uniqueEdgeByElem(:,iEdge);
                eInElem(:,iEdge) = obj.uniqueToAll(edge);
            end
            obj.edgesInElem = eInElem;
        end

        function edgeByElem = computeUniqueEdgesByElem(obj)
            edges(:,1) = 1:obj.nAllEdges;
            edgeByElem = reshape(edges,obj.nEdgeByElem,obj.nElem);
            edgeByElem = edgeByElem';
        end

        function computeLocalOrientedEdgeConnec(obj)
            edgeConnec = zeros(obj.nElem,obj.nEdgeByElem,obj.nNodeByEdge);
            for iedge = 1:obj.nEdgeByElem
                isOriented = obj.isEdgeOriented(iedge);
                nodeA = obj.localEdgesInElem(iedge,1);
                nodeB = obj.localEdgesInElem(iedge,2);
                edgeConnec(isOriented,iedge,1)  = nodeA;
                edgeConnec(isOriented,iedge,2)  = nodeB;
                edgeConnec(~isOriented,iedge,1) = nodeB;
                edgeConnec(~isOriented,iedge,2) = nodeA;
            end
            obj.localNodeByEdgeByElem = edgeConnec;
        end

        function itIs = isEdgeOriented(obj,iedge)
            globalNodes = obj.computeGlobalNode(iedge);
            localNodes  = obj.computeLocalNodes(iedge);
            itIs = obj.isEquallyOriented(globalNodes,localNodes);
        end

        function globalNodes = computeGlobalNode(obj,iedge)
            edges       = obj.edgesInElem(:,iedge);
            globalNodes = obj.nodesInEdges(edges,:);
        end

        function localNodes = computeLocalNodes(obj,iedge)
            localNode   = obj.localEdgesInElem(iedge,:);
            localNodes  = obj.nodesByElem(:,localNode);
        end

        function edges = computeAllEdgesOfVertex(obj,vertex)
            vertexInEdges = obj.nodesInEdges;
            isInEdge = vertexInEdges == vertex;
            allEdges  = 1:size(isInEdge,1);
            allEdgesA(:,1) = allEdges(isInEdge(:,1));
            allEdgesB(:,1) = allEdges(isInEdge(:,2));
            edges = [allEdgesA;allEdgesB];
        end

        function oV = computeOtherVertexOfEdge(obj,edge,vertex)
            vertexesInEdges = obj.nodesInEdges;
            vertexesOfEdge  = vertexesInEdges(edge,:);
            oVB = setdiff(vertexesOfEdge(:,2),vertex);
            oVA = setdiff(vertexesOfEdge(:,1),vertex);
            oV = [oVB;oVA];
        end

    end

    methods (Access = private, Static)

        function itIs = isEquallyOriented(gNode,lNode)
            itIs = sum(abs(gNode - lNode),2) < 1;
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [cl] = qt2cl(qt)
% Convert volume from US liquid quarts to centiliters.
% Chad Greene 2012
cl = qt*94.6352946;


text:
%% RATE OF CONVERGENCE OF MIXED FINITE ELEMENT METHOD (RT0-P0) FOR POISSON EQUATION
%
% This example is to show the rate of convergence of RT0-P0 mixed finite
% element approximation of the Poisson equation on the unit square with the
% following boundary conditions:
%
% - pure Dirichlet boundary condition.
% - Pure Neumann boundary condition.
% - mxied boundary condition.
%
% The basis, data structure and numerical test is summarized in <a
% href="matlab:ifem Poisson3RT0femrate">Poisson3RT0femrate</a>.
%
% The optimal rates of convergence for u and sigma are observed, namely,
% 1st order for L2 norm of u, L2 norm of sigma and H(div) norm of sigma.
% The 2nd order convergent rates between two discrete functions ||uI-uh||
% and ||sigmaI-sigmah|| are known as superconvergence.
%
% Triangular preconditioned GMRES (default solver) and Uzawa preconditioned
% CG converges uniformly in all cases. Traingular preconditioner is two
% times faster than PCG although GMRES is used.
% See also PoissonBDM1mfemrate, Poissonfemrate, Poissonafemrate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Setting
[node,elem] = cubemesh([0,1,0,1,0,1],0.5);
mesh = struct('node',node,'elem',elem);
option.L0 = 1;
option.maxIt = 3;
option.printlevel = 1;
option.elemType = 'RT0';
% pde = mixedPossiondata;

%% Dirichelt for u and Neumann boundary condition for sigma
option.solver = 'uzawapcg';
pde = sincosdata3;
mesh.bdFlag = setboundary3(node,elem,'Dirichlet');
mfemPoisson3(mesh,pde,option);

%% Neumann for u and Dirichlet boundary condition for sigma.
option.solver = 'tri';
pde = sincosdata3;
mesh.bdFlag = setboundary3(node,elem,'Neumann');
mfemPoisson3(mesh,pde,option);

repo: ifem


text:
function X = Yi2X(Y, i)
if i == 3
    X = shiftdim(Trans_Faces(Y),i+1);
elseif i == 2
    X = shiftdim(Trans_Faces(Y),i);
else
    X = shiftdim(Trans_Faces(Y),i-1);
end
end

author: andrewssobral


text:
function r = sqrt(a)
%SQRT         Hessian (elementwise) square root
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

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

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = sqrt(a.x);
    rx1 = 0.5 / r.x;
    r.dx = a.dx * rx1;
    r.hx = ( a.hx - reshape( ((0.25/a.x)*a.dx) * a.dx.',size(a.hx)) ) * rx1;

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = sqrt(a.x);
    if issparse(a.hx)               % input sparse

      ax = 0.5./full(r.x(:));
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        rdx = ax(ja).*sa(:);
        adx1 = (-0.25) ./ a.x(ja);
        if isa(a.x,'intval')          % sparse intval
          adx1 = times(adx1(:),rdx,0);
          if isreal(rdx)
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          else
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          end
          if adx1.complex
            adx1 = intval( sparse(ia,ja,adx1.mid,N,sizeax) , sparse(ia,ja,adx1.rad,N,sizeax) , 'midrad' );
          else
            adx1 = intval( sparse(ia,ja,adx1.inf,N,sizeax) , sparse(ia,ja,adx1.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,rdx,N,sizeax);
          adx1 = sparse(ia,ja,adx1(:).*rdx,N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,adx1,a.dx);
      end
      [ia,ja,sa] = find(a.hx);      % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      r.x = sqrt(a.x);
      rx1 = 0.5 ./ ( r.x(:).' );
      rx1 = rx1(ones(N*N,1),:);
      r.dx = a.dx .* rx1(1:N,:);
      adx = repmat(0.25./(a.x(:).'),N,1) .* a.dx;
      r.hx = ( a.hx - adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:) ) .* rx1;

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
% GetSegyTraceData : Get Segy trace data if filehandle
%
% Call :
%
%   tracedata=GetSegyTraceData(segyid,ns,SegyHeader)
%

%
% (C) 2001-2004 Thomas Mejer Hansen, tmh@gfy.ku.dk/thomas@cultpenguin.com
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
function tracedata=GetSegyTraceData(segyid,ns,SegyHeader,SkipData,it)
%
% Get Segy Trace Data.
%
%

SkipData=[];

if nargin==0,
    SegymatVerbose(exist('segyid'))
    SegymatVerbose([mfilename,' : SEGYID not specified - exiting'])
    tracedata=[];
    return
end
if nargin==1
  SegymatVerbose([mfilename,' : NS not specified - exiting'])
  tracedata=[];
  return
end
if nargin==2
    SegyHeader.DataFormat='float32';
    SegyHeader.BytesPerSample=32;
    SegyHeader.DataSampleFormat=5; % IEEE
    SegymatVerbose(['Dataformat not specified -  dataformat->',SegyHeader.DataFormat])
end;

if isempty(SkipData)==1, SkipData=0; end

Revision=SegyHeader.SegyFormatRevisionNumber;
if Revision>0, Revision=1; end
Format=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).format;

BPS=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).bps;
if (SkipData==1)
    SkipBytes=ns*BPS/8;
    fseek(segyid,SkipBytes,'cof');
    tracedata=[];
else
    % SegymatVerbose([mfilename,' : ',Format,'. ns=',num2str(ns)])
    try

       if nargin>4

           % Skip to the it'th sample
           SkipBytes=(it-1)*BPS/8;
           fseek(segyid,SkipBytes,'cof');

           tracedata=fread(segyid,1,Format);

           SkipBytes=(ns-it)*BPS/8;
           fseek(segyid,SkipBytes,'cof');

       else
         tracedata=fread(segyid,ns,Format);
       end
    catch
      SegymatVerbose([mfilename,' : Error using fread - Possibly ''ns'' is negative -' ...
	    ' check byteorder-'])
      tracedata=[];
    end

    if (strcmp(Format,'uint32')==1), % IBM FLOATING POINT
        % CONVERT FROM FLOATING POINT
        verbose=1;
        if verbose>1, SegymatVerbose([mfilename,'Converting from IBM, DataFormat :',SegyHeader.DataFormat]); end
        try
          tracedata=ibm2num(uint32(tracedata));
        catch
          % SegymatVerbose([mfilename,' : SOMETHING BAD HAPPENED WHEN CONVERTING FROM IBM FLOATS TO IEEE. ARE YOU SURE DATA ARE IBM FLOAT FORMATTED ?' ])
          % tracedata=0.*tracedata;
          % return

        end
    end;

end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function results = YahooGeoCodeAddress(appId, location)
% YAHOOGEOCODEADDRESS given a developer Application ID (appId)
% and string representing an address LOCATION, returns
% a structure containing the possible matches for the location. This
% function uses the Yahoo Maps Service - Geocoding API
% For more information please visit to get an appId
% http://developer.yahoo.com/maps/rest/V1/geocode.html

% Examples
% This is only a demo AppId 'YD-9G7bey8_JXxQP6rxl.fBFGgCdNjoDMACQA--'
% invalid data may result in using this AppId, you must get a real appId.
%
% results = GeoCodeAddress('YD-9G7bey8_JXxQP6rxl.fBFGgCdNjoDMACQA--','The MathWorks, Natick MA')
% results = GeoCodeAddress('YD-9G7bey8_JXxQP6rxl.fBFGgCdNjoDMACQA--','3 Apple Hill drive, Natick MA')
% results = GeoCodeAddress('YD-9G7bey8_JXxQP6rxl.fBFGgCdNjoDMACQA--','Natick MA')

% Throw a warning if they are using the demo AppId
if strcmp(appId,'YD-9G7bey8_JXxQP6rxl.fBFGgCdNjoDMACQA--')
   warning('YAHOOGEOCODE:usingdemoappid',...
       ['You are using the demo AppId provided by yahoo.', char(13),...
       'To get a valid AppID you must register at',char(13),...
       'http://developer.yahoo.com/maps/rest/V1/geocode.html']);
end

% The URL for Yahoos GeoCode API
url = 'http://local.yahooapis.com/MapsService/V1/geocode';

% Replace spaces with + signs, so it can go inside a URL
location = regexprep(location,'\s+','+');

% Query yahoo for the data
try
    docNode = xmlread([url '?appid=' appId '&location=' location]);
catch  %#ok<CTCH>
    error('YAHOOGEOCODE:urlreaderror',...
        ['Could not reach Yahoo server, please check your AppID ', char(13),...
        'to make sure it is valid.']);
end

% Check for errors
if docNode.getElementsByTagName('Error').getLength>0
   % Display the first error
   errorMessage = ...
       char(docNode.getElementsByTagName('Error').item(0).getTextContent);
   error('YAHOOGEOCODE:yahooerror',errorMessage);
end

% Loop over the results and create a struct for them
resultList = docNode.getElementsByTagName('Result');
results = [];
for idx = 1:resultList.getLength
    resultNode = resultList.item(idx-1);
    results(idx).precision = char(resultNode.getAttribute('precision')); %#ok<*AGROW>
    results(idx).Latitude = GetElementText(resultNode,'Latitude');
    results(idx).Longitude = GetElementText(resultNode,'Longitude');
    results(idx).Address = GetElementText(resultNode,'Address');
    results(idx).City = GetElementText(resultNode,'City');
    results(idx).State = GetElementText(resultNode,'State');
    results(idx).Zip = GetElementText(resultNode,'Zip');
    results(idx).Country = GetElementText(resultNode,'Country');
end
end

function elementText = GetElementText(resultNode,elementName)
% GETELEMENTTEXT given a result node and an element name
% returns the text within that node

elementText = ...
    resultNode.getElementsByTagName(elementName).item(0).getTextContent;

% remove any rogue commas and casts to a char
elementText = regexprep(char(elementText),',','');
end


function gradx = singletrack_gradx_f(in1,in2,in3)
%SINGLETRACK_GRADX_F
%    GRADX = SINGLETRACK_GRADX_F(IN1,IN2,IN3)

%    This function was generated by the Symbolic Math Toolbox version 8.3.
%    09-Jan-2020 11:52:33

I_z = in3(2,:);
M = in3(1,:);
T = in2(2,:);
V_vx = in1(4,:);
V_vy = in1(5,:);
c_f = in3(8,:);
c_r = in3(9,:);
delta = in2(1,:);
deltamax = in3(5,:);
l_f = in3(3,:);
l_r = in3(4,:);
maxbrakeWForce = in3(6,:);
psi_dot = in1(6,:);
vpsi = in1(3,:);
t2 = cos(vpsi);
t3 = sin(vpsi);
t4 = l_f.*psi_dot;
t5 = l_r.*psi_dot;
t6 = V_vx.^2;
t7 = 1.0./I_z;
t8 = 1.0./M;
t9 = T.*5.0e+1;
t10 = delta.*5.0e+1;
t11 = deltamax.*5.0e+1;
t14 = V_vx.*1.0e+2;
t12 = -t5;
t13 = -t9;
t15 = -t10;
t16 = -t11;
t17 = V_vy+t4;
t18 = tanh(t14);
t19 = t17.^2;
t20 = V_vy+t12;
t21 = t13+5.0e+1;
t22 = t18.^2;
t23 = t13-5.0e+1;
t27 = t11+t15;
t31 = t15+t16;
t24 = t20.^2;
t25 = exp(t21);
t26 = exp(t23);
t28 = exp(t27);
t29 = t6+t19;
t32 = exp(t31);
t34 = t22.*1.0e+2;
t30 = t25+1.0;
t33 = t26+1.0;
t35 = t6+t24;
t36 = t28+1.0;
t38 = t32+1.0;
t39 = 1.0./t29;
t44 = t34-1.0e+2;
t37 = 1.0./t30;
t40 = 1.0./t33;
t41 = 1.0./t36;
t42 = 1.0./t35;
t43 = 1.0./t38;
t45 = t37.*5.0e+1;
t46 = t37-1.0;
t47 = t40.*5.0e+1;
t49 = deltamax.*t41;
t50 = t41-1.0;
t52 = V_vx.*c_r.*l_r.*t42;
t53 = t43-1.0;
t48 = -t45;
t51 = -t47;
t54 = deltamax.*t53;
t55 = T.*t40.*t46;
t57 = t9.*t40.*t46;
t58 = delta.*t43.*t50;
t56 = -t55;
t59 = -t58;
t61 = t48+t51+t57+5.0e+1;
t60 = t40+t46+t56;
t62 = exp(t61);
t64 = t49+t54+t59;
t63 = t62+1.0;
t65 = sin(t64);
t67 = cos(t64);
t66 = 1.0./t63;
t69 = V_vx.*c_f.*l_f.*t39.*t67;
t68 = t66-1.0;
gradx = reshape([0.0,0.0,-V_vx.*t3-V_vy.*t2,t2,-t3,0.0,0.0,0.0,0.0,V_vx.*t2-V_vy.*t3,t3,t2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,t8.*(-c_f.*t17.*t39.*t65+(maxbrakeWForce.*t44.*t60.*t68)./2.0+(maxbrakeWForce.*t44.*t60.*t67.*t68)./2.0),t8.*(psi_dot+V_vx.*c_f.*t39.*t65),t8.*(V_vy+V_vx.*c_f.*l_f.*t39.*t65),0.0,0.0,0.0,0.0,t8.*(c_r.*t20.*t42+c_f.*t17.*t39.*t67+(maxbrakeWForce.*t44.*t60.*t65.*t68)./2.0),-t8.*(psi_dot+V_vx.*c_r.*t42+V_vx.*c_f.*t39.*t67),-t8.*(V_vy-t52+t69),0.0,0.0,0.0,0.0,t7.*(-c_r.*l_r.*t20.*t42+c_f.*l_f.*t17.*t39.*t67+(l_f.*maxbrakeWForce.*t44.*t60.*t65.*t68)./2.0),t7.*(t52-t69),-t7.*(l_f.*t69+l_r.*t52),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[7,7]);


text:
function [isih] = ft_spike_isi(cfg, spike)

% FT_SPIKE_ISI computes the interspike interval histogram
%
% The input SPIKE should be organised as
% a) the spike datatype, obtained from FT_SPIKE_MAKETRIALS
% b) the raw datatype, containing binary spike trains, obtained from
% FT_APPENDSPIKE or FT_CHECKDATA. In this case the raw datatype is
% converted to the spike datatype.
%
% Use as
%   [isih] = ft_spike_isi(cfg, spike)
%
% Configurations:
%   cfg.outputunit       = 'spikecount' (default) or 'proportion' (sum of all bins = 1).
%   cfg.spikechannel     = string or index of spike channels to
%                          trigger on (default = 'all')
%                          See FT_CHANNELSELECTION for details.
%   cfg.trials           = numeric selection of trials (default = 'all')
%   cfg.bins             = ascending vector of isi bin edges.
%   cfg.latency          = [begin end] in seconds, 'max' (default), 'min', 'prestim'(t<=0), or
%                          'poststim' (t>=0).
%                          If 'max', we use all available latencies.
%                          If 'min', we use only the time window contained by all trials.
%                          If 'prestim' or 'poststim', we use time to or
%                          from 0, respectively.
%`  cfg.keeptrials       = 'yes' or 'no'. If 'yes', we keep the individual
%                           isis between spikes and output as isih.isi
%   cfg.param            = string, one of
%      'gamfit'      : returns [shape scale] for gamma distribution fit
%      'coeffvar'    : coefficient of variation (sd / mean)
%      'lv'          : Shinomoto's Local Variation measure (2009)
%
% Outputs:
%   isih.avg             = nUnits-by-nBins interspike interval histogram
%   isih.time            = 1 x nBins bincenters corresponding to isih.avg
%   isih.isi             = 1-by-nUnits cell with interval to previous spike per spike.
%                          For example isih.isi{1}(2) = 0.1 means that the
%                          second spike fired was 0.1 s later than the
%                          first. Note that jumps within trials or first
%                          spikes within trials are given NaNs.
%   isih.label           = 1-by-nUnits cell-array with labels

% Copyright (C) 2010-2012, Martin Vinck
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance spike

% check if data is of proper format
spike = ft_checkdata(spike,'datatype', 'spike', 'feedback', 'yes');

% get the default options
cfg.outputunit   = ft_getopt(cfg,'outputunit','spikecount');
cfg.spikechannel = ft_getopt(cfg,'spikechannel', 'all');
cfg.trials       = ft_getopt(cfg,'trials', 'all');
cfg.latency      = ft_getopt(cfg,'latency','maxperiod');
cfg.keeptrials   = ft_getopt(cfg,'keeptrials', 'yes');
cfg.bins         = ft_getopt(cfg,'bins', 0:0.002:1);
cfg.param        = ft_getopt(cfg,'param', 'coeffvar');

% ensure that the options are valid
cfg = ft_checkopt(cfg,'outputunit','char', {'spikecount', 'proportion'});
cfg = ft_checkopt(cfg,'bins', 'ascendingdoublevector');
cfg = ft_checkopt(cfg,'spikechannel',{'cell', 'char', 'double'});
cfg = ft_checkopt(cfg,'latency', {'char', 'ascendingdoublebivector'});
cfg = ft_checkopt(cfg,'trials', {'char', 'doublevector', 'logical'});
cfg = ft_checkopt(cfg,'keeptrials', 'char', {'yes', 'no'});
cfg = ft_checkopt(cfg,'param', 'char', {'gamfit', 'coeffvar', 'lv'});

cfg = ft_checkconfig(cfg, 'allowed', {'param', 'outputunit', 'bins', 'spikechannel', 'latency', 'trials', 'keeptrials'});

% get the number of trials or change DATA according to cfg.trials
if  strcmp(cfg.trials,'all')
  cfg.trials = 1:size(spike.trialtime,1);
elseif islogical(cfg.trials) || all(cfg.trials==0 | cfg.trials==1)
  cfg.trials = find(cfg.trials);
end
cfg.trials = sort(cfg.trials);

cfg.channel = ft_channelselection(cfg.spikechannel, spike.label);
spikesel    = match_str(spike.label, cfg.channel);
nUnits      = length(spikesel); % number of spike channels
if nUnits==0, error('No spikechannel selected by means of cfg.spikechannel'); end

% determine the duration of each trial
begTrialLatency = spike.trialtime(cfg.trials,1);
endTrialLatency = spike.trialtime(cfg.trials,2);

% select the latencies
if strcmp(cfg.latency,'minperiod')
  cfg.latency = [max(begTriallatency) min(endTriallatency)];
elseif strcmp(cfg.latency,'maxperiod')
  cfg.latency = [min(begTrialLatency) max(endTrialLatency)];
elseif strcmp(cfg.latency,'prestim')
  cfg.latency = [min(begTrialLatency) 0];
elseif strcmp(cfg.latency,'poststim')
  cfg.latency = [0 max(endTrialLatency)];
end

% construct the isi bins
bins  = cfg.bins;
nBins = length(bins)-1;

% compute the interspike interval
keepTrials = strcmp(cfg.keeptrials,'yes');
if keepTrials, isiSpike = cell(1,nUnits); end  % contains the individual histogram spike times
isihist  = zeros(nUnits,nBins+1); % isi histogram
out = [];
for iUnit = 1:nUnits
  unitIndx = spikesel(iUnit);

  % only select the spikes in the right latencies and trials
  spikeTrials    = spike.trial{unitIndx}(:)'; % ensure row vector
  spikeTimes     = spike.time{unitIndx}(:)'; % ensure row vector

  % select only the spikes in the window and with the selected trials
  spikesInTrials = ismember(spikeTrials, cfg.trials); % row vec
  spikesInWin    = spikeTimes>=cfg.latency(1)&spikeTimes<=cfg.latency(2); % row vec
  spikeTimes     = spikeTimes(spikesInTrials & spikesInWin);
  spikeTrials    = spikeTrials(spikesInTrials & spikesInWin);

  % find the spikes that jumped to the next trial, replace them with NaNs
  trialJump = logical([1 diff(spikeTrials(:)')]);
  isi = [NaN diff(spikeTimes(:)')];
  isi(trialJump) = NaN;

  switch cfg.param
  case 'coeffvar'
    out(iUnit) = nanstd(isi)./nanmean(isi);
  case 'gamfit'
    data = isi(~isnan(isi));  % remove the nans from isiSpike
    if ~isempty(data)
      [out(iUnit,:)] = mle(data,'distribution', 'gamma'); % fit a gamma distribution
    else
      out(iUnit,:)   = [NaN NaN];
    end
  case 'lv'
    dIsi     = isi(1:end-1) - isi(2:end);
    sumIsi   = isi(1:end-1) + isi(2:end);
    sl       = ~isnan(dIsi) & ~isnan(sumIsi); % remove if one has nan
    df       = sum(sl) - 1;
    out(iUnit) = (3/df)*sum((dIsi(sl)./sumIsi(sl)).^2);
  end

  % convert and store the isi
  if keepTrials, isiSpike{iUnit} = isi; end
  isihist(iUnit,:)   = histc(isi,bins);
end

isihist(:,end) = []; % the last number is only an equality to a bin edge
if strcmp(cfg.outputunit,'proportion')
  isihist = isihist./repmat(nansum(isihist,2),1,size(isihist,2));
end

% gather the rest of the results
if strcmp(cfg.keeptrials,'yes'), isih.isi  = isiSpike; end
isih.time        = bins(1:end-1);
isih.avg         = isihist;
isih.dimord      = 'chan_time';
isih.label       = spike.label(spikesel);
param = cfg.param;
isih.(param) = out;

% do the general cleanup and bookkeeping at the end of the function

ft_postamble previous   spike
ft_postamble provenance isih
ft_postamble history    isih

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% FUNCTION split_data
%   Splitting multi-task data into training / testing by percentage.
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   percent: percentage of the splitting range (0, 1)
%
%% OUTPUT
%   X_train: the split of X that has the specified percent of samples
%   Y_train: the split of Y that has the specified percent of samples
%   X_test: the split of X that has the remaining samples
%   Y_test: the split of Y that has the remaining samples
%   selIdx: the selection index of for X_train and Y_train for each task
%%

function [X_train, Y_train, X_test, Y_test, selIdx] = split_data(X, Y, percent)

if percent > 1 || percent < 0
    error('splitting percentage error')
end

task_num = length(X);

selIdx = cell(task_num, 0);
X_train = cell(task_num, 0);
Y_train = cell(task_num, 0);
X_test = cell(task_num, 0);
Y_test = cell(task_num, 0);

for t = 1:task_num
    task_sample_size = length(Y{t});
    tSelIdx = randperm(task_sample_size) < task_sample_size * percent;

    selIdx{t} = tSelIdx;

    X_train{t} = X{t}(tSelIdx,:);
    Y_train{t} = Y{t}(tSelIdx,:);
    X_test{t} = X{t}(~tSelIdx,:);
    Y_test{t} = Y{t}(~tSelIdx,:);

end

author: gingsmith

repo: fmtl

save_path: github-repos/MATLAB/gingsmith-fmtl


text:
Vh = zeros(size(mesh.t,1),3);

feEvalBas = @EvalNed1Bas3D;
IfeEvalBas = @evalNed1IFEBas3D;
option.ng = 1;
option.nge = 1;

p = mesh.p; t = mesh.t;
fem = genNedFEM3D(mesh,bc,option);
ntID = find(mesh.tLoc>0);
gxN = fem.gx(ntID,:); gyN = fem.gy(ntID,:); gzN = fem.gz(ntID,:);

dof = fem.ldof;
ntID = find(mesh.tLoc>0);
uhK = uh(fem.g2ldof);
uK1 = 0; uK2 = 0; uK3 = 0;
for i = 1:dof
    uK1 = uK1 + uhK(ntID,i).*feEvalBas(fem.bas, ntID, gxN, gyN, gzN, i, 0, 1).*fem.t_e_orit(ntID,i);
    uK2 = uK2 + uhK(ntID,i).*feEvalBas(fem.bas, ntID, gxN, gyN, gzN, i, 0, 2).*fem.t_e_orit(ntID,i);
    uK3 = uK3 + uhK(ntID,i).*feEvalBas(fem.bas, ntID, gxN, gyN, gzN, i, 0, 3).*fem.t_e_orit(ntID,i);
end

Vh(ntID,:) = [uK1,uK2,uK3];

% interface elements
epsm = 8.85*10^(-3)*5;
epsp = 8.85*10^(-3);
sigm = 100;
sigp = 1;
mum = (4*pi);
mup = (4*pi)*2;
bm = epsm/deltaT^2 + sigm/(2*deltaT);
bp = epsp/deltaT^2 + sigp/(2*deltaT);
am = mum^(-1); ap = mup^(-1);
femI = genNed1IFEM3D(mesh,fem,am,ap,bm,bp,option);
gx = femI.gx; gy = femI.gy; gz = femI.gz;
tIntfID = femI.tIntfID;
itID = find(mesh.tLoc < 0); ntI = length(itID);
uhK = uh(fem.g2ldof); uhKitmp = uhK(itID,:); uhKi = zeros(length(femI.g2ldof),6);
uK1 = 0; uK2 = 0; uK3 = 0;
id = 0;
for i = 1:ntI
    tmp = uhKitmp(i,:);
    uhKi(id+1:id+femI.quadT(i),:) = repmat(tmp,femI.quadT(i),1);
    id = id + femI.quadT(i);
end

for i = 1:dof
    BAS = femI.bas(:,:,i);
    uK1 = uK1 + uhKi(:,i).*IfeEvalBas(BAS, gx, gy, gz, 0, 1).*fem.t_e_orit(tIntfID,i);
    uK2 = uK2 + uhKi(:,i).*IfeEvalBas(BAS, gx, gy, gz, 0, 2).*fem.t_e_orit(tIntfID,i);
    uK3 = uK3 + uhKi(:,i).*IfeEvalBas(BAS, gx, gy, gz, 0, 3).*fem.t_e_orit(tIntfID,i);
end
tIntfIDtmp = -mesh.tLoc(tIntfID);
Vi1 = accumarray(tIntfIDtmp,uK1)./femI.quadT;
Vi2 = accumarray(tIntfIDtmp,uK2)./femI.quadT;
Vi3 = accumarray(tIntfIDtmp,uK3)./femI.quadT;
Vh(itID,:) = [Vi1,Vi2,Vi3];

fMall = (mesh.p(mesh.t(:,1),:)+mesh.p(mesh.t(:,2),:)+...
    mesh.p(mesh.t(:,3),:)+mesh.p(mesh.t(:,4),:))/4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol = 10^(-12);
zCut = -0.3125;
fp1 = mesh.p(mesh.f(:,1),:);
fp2 = mesh.p(mesh.f(:,2),:);
fp3 = mesh.p(mesh.f(:,3),:);
fZid = (abs(fp1(:,3)-zCut)<tol).*(abs(fp2(:,3)-zCut)<tol).*(abs(fp3(:,3)-zCut)<tol);
fZid = find(fZid==1);
fmpt = (fp1(fZid,:)+fp2(fZid,:)+fp3(fZid,:))/3;
fmptPlan = fmpt(:,1:2);
zDT = delaunay(fmptPlan);

fZTid = mesh.f_t(fZid,1);
VhCut = Vh(fZTid,:);

%%%%% trisulf plot
trisurf(zDT,fmptPlan(:,1),fmptPlan(:,2),VhCut(:,1));
shading interp

%%%%%% contour plot
tricontour(fmptPlan,zDT,VhCut(:,1),10)

%%%%% vector field plot
quiver3(fmpt(:,1),fmpt(:,2),fmpt(:,3),VhCut(:,1),VhCut(:,2),VhCut(:,3))

tid1 = find(mesh.tLoc==1);
tid1 = tid1(1:50:end);
quiver3(fMall(tid1,1),fMall(tid1,2),fMall(tid1,3),Vh(tid1,1),Vh(tid1,2),Vh(tid1,3),20)

tid1 = find(mesh.tLoc==1);
tid1 = tid1(1:50:end);
quiver5(fMall(tid1,1),fMall(tid1,2),fMall(tid1,3),...
    Vh(tid1,1),Vh(tid1,2),Vh(tid1,3),5);

tidzn = find(fMall(:,3)<0);
tidzn = tidzn(1:100:end);
quiver5(fMall(tidzn,1),fMall(tidzn,2),fMall(tidzn,3),...
    Vh(tidzn,1),Vh(tidzn,2),Vh(tidzn,3),3);


text:
function range = p00_range ( prob )

%*****************************************************************************80
%
%% P00_RANGE returns an interval bounding the root for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the number of the problem.
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  if ( prob == 1 )
    range = p01_range ( );
  elseif ( prob == 2 )
    range = p02_range ( );
  elseif ( prob == 3 )
    range = p03_range ( );
  elseif ( prob == 4 )
    range = p04_range ( );
  elseif ( prob == 5 )
    range = p05_range ( );
  elseif ( prob == 6 )
    range = p06_range ( );
  elseif ( prob == 7 )
    range = p07_range ( );
  elseif ( prob == 8 )
    range = p08_range ( );
  elseif ( prob == 9 )
    range = p09_range ( );
  elseif ( prob == 10 )
    range = p10_range ( );
  elseif ( prob == 11 )
    range = p11_range ( );
  elseif ( prob == 12 )
    range = p12_range ( );
  elseif ( prob == 13 )
    range = p13_range ( );
  elseif ( prob == 14 )
    range = p14_range ( );
  elseif ( prob == 15 )
    range = p15_range ( );
  elseif ( prob == 16 )
    range = p16_range ( );
  elseif ( prob == 17 )
    range = p17_range ( );
  elseif ( prob == 18 )
    range = p18_range ( );
  elseif ( prob == 19 )
    range = p19_range ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_RANGE - Fatal error%\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', prob );
    error ( 'P00_RANGE - Fatal error%' );
  end

  return
end


text:
function xy_io_test01 ( )

%*****************************************************************************80
%
%% TEST01 tests XY_EXAMPLE, XY_WRITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2009
%
%  Author:
%
%    John Burkardt
%
  point_num = 300;
  filename = 'xy_io_test01.xy';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  XY_EXAMPLE sets up sample XY data.\n' );
  fprintf ( 1, '  XY_WRITE writes an XY file.\n' );

  xy = xy_example ( point_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XY_EXAMPLE has created the data.\n' );

  xy_write ( filename, point_num, xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XY_WRITE wrote the header and data for %s\n', filename );
  fprintf ( 1, '  Number of points = %d\n', point_num );

  return
end


text:
classdef TestGetOptimalDFTSize
    %TestGetOptimalDFTSize

    methods (Static)
        function test_1
            N = cv.getOptimalDFTSize(101);
            validateattributes(N, {'numeric'}, {'scalar', 'integer'});
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function gX = rbfperiodicKernDiagGradX(kern, X)

% RBFPERIODICKERNDIAGGRADX Gradient of RBFPERIODIC kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the RBF derived periodic kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : rbfperiodicKernParamInit, kernDiagGradX, rbfperiodickernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

gX = zeros(size(X));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef ThrottleTermCondOptimVar < AbstractOptimizationVariable
    %ThrottleTermCondOptimVar Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) ThrottleTermCondition = ThrottleTermCondition(0);

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = ThrottleTermCondOptimVar(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.targetThrottle;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(obj.useTf);
            ub = obj.ub(obj.useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = obj.useTf;
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = useTf;
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.targetThrottle = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            nameStrs = {sprintf('Event %i Throttle Termination Condition', evtNum)};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ objects, nElems ] = buildGroundTruth( objects, threshold_detection, rebuild )
%BUILDGROUNDTRUTH Analyzes the GT and applies a matching on each of the
% object candidate windows.

    nSamples = length(objects);
    nElems = 0;

    for j = 1:nSamples
        nGT = length(objects(j).ground_truth);

        if(rebuild)

            %% Initialize object candidates
            nObjs = length(objects(j).objects);
            for k = 1:nObjs
                objects(j).objects(k).OS = zeros(1, nGT);
                objects(j).objects(k).trueLabel = 'No Object';
                objects(j).objects(k).trueLabelId = [];
            end

            %% For each object
            for k = 1:nGT
                %% Check if any object in objects(j).objects matches
                %  with the ground_truth for assign them the true label!
                GT = objects(j).ground_truth(k);
                GT.height = (GT.BRy - GT.ULy + 1);
                GT.width = (GT.BRx - GT.ULx + 1);
                GT.area = GT.height * GT.width;

                %% Check for each object candidate, if it fits the current true object
                count_candidate = 1;
                for w = objects(j).objects
                    if(~isempty(w.ULx)) % if the list of candidates is not empty
                        % Check area and intersection on current window "w"
                        w.height = (w.BRy - w.ULy + 1);
                        w.width = (w.BRx - w.ULx + 1);
                        w.area = w.height * w.width;

                        % Check intersection
    %                     count_intersect_old = rectint([GT.ULy, GT.ULx, GT.height, GT.width], [w.ULy, w.ULx, w.height, w.width]);
                        x_overlap = max(0, min(GT.BRx,w.BRx) - max(GT.ULx,w.ULx));
                        y_overlap = max(0, min(GT.BRy,w.BRy) - max(GT.ULy,w.ULy));
                        count_intersect = x_overlap * y_overlap;

                        % Calculate overlap score
                        OS = count_intersect / (GT.area + w.area - count_intersect);

                        if(OS >= threshold_detection) % object detected!
                            label = objects(j).ground_truth(k).name;
                            % If OS bigger than previous, then assign this
                            if(max(w.OS) < OS)
                                w.trueLabel = label;
                                w.trueLabelId = k;
                            end
                        end
                        w.OS(k) = OS;

                        % Store w object
                        objects(j).objects(count_candidate).OS = w.OS;
                        objects(j).objects(count_candidate).trueLabel = w.trueLabel;
                        objects(j).objects(count_candidate).trueLabelId = w.trueLabelId;

                        count_candidate = count_candidate + 1;
                    end
                end
                nElems = nElems+count_candidate-1;
            end
        else % Only pick the maximum OS (already calculated) if it is over the given threshold
            count_candidate = 1;
            for w = objects(j).objects
                if(w.OS(w.trueLabelId) < threshold_detection)
                    objects(j).objects(count_candidate).trueLabel = 'No Object';
                    objects(j).objects(count_candidate).trueLabelId = [];
                end
                count_candidate = count_candidate + 1;
            end
            nElems = nElems+count_candidate-1;
        end
    end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function F = initBPHMMSamplerFromScratch_FeatureMatrix( data_struct, initParams,outParams)
%  Features can be initialized in a few ways
%    1) unique set for each obj
%    2) unique set for each category
%    3) shared set for all objects
%           either from IBP prior
%                    or given specific number K of features to use

nObj = length( data_struct );

if isfield( initParams.F, 'nTotal' )
    K =  initParams.F.nTotal;
    F = ones( nObj, K );

    if outParams.doPrintHeaderInfo
        fprintf( '\t F : %d global features shared by all objects \n', K );
    end

elseif isfield( initParams.F, 'nUniquePerObj' )
    Kii = initParams.F.nUniquePerObj;
    F = zeros( nObj, Kii * nObj );

    for ii = 1:nObj
        F(ii, (ii-1)*Kii + [1:Kii] ) = 1;
    end
    if outParams.doPrintHeaderInfo
        fprintf( '\t F : each obj. assigned %d unique features \n', Kii );
    end

end

author: michaelchughes

repo: NPBayesHMM


text:
function googlesearch_test
%GOOGLESEARCH_TEST Tests for GOOGLESEARCH.

% Make a simple call.
s = googlesearch('matlab');
if ~isstruct(s)
    error('Did not return a structure.')
end
if ~isfield(s,'ResultElements')
    error('No field "ResultElements".')
end
resultElements = s.ResultElements;
if ~isstruct(resultElements)
    error('ResultElements is not a structure."')
end
if ~isequal(size(resultElements),[1 10])
    error('"ResultElements is not the right size.')
end
url = resultElements(5).URL;
if ~ischar(url) & ~isempty(url)
    error('URLs not returned properly.')
end

%% Make a fancier call.
s = googlesearch('matlab mathworks',90);
%%
if ~isequal(size(s.ResultElements),[1 90])
    error('"ResultElements is not the right size.')
end
if (s.StartIndex ~= 1)
    error('StartIndex is wrong.')
end
if (s.EndIndex ~= 90)
    error('EndIndex is wrong.')
end
title = s.ResultElements(72).Title;
if ~ischar(title) && ~isempty(title)
    error('"Title" is wrong.')
end


text:
function [rmses, av_rmse] = rmse_f(x,y)

% RMSE calculator

% 02/12/2009 Version 1.5    - Hyperspectral support
% 05/03/2010 Version 2.0    - Double precision, RAM efficiency
% 25/06/2010 Version 2.2    - Excel Output option
% 08/08/2011 Version 2.4F   - Function Version

% Author: Aristidis D. Vaiopoulos

sizi = size(x);
if max(size(size(x))) == 2
    bands = 1;
else
    bands = sizi(1,3);
end

nres = sizi(1,1)*sizi(1,2);

rmses = zeros(bands,1);
for i = 1:bands
    xt = double(x(:,:,i));
    yt = double(y(:,:,i));
    rmses(i) = sqrt( (sum(sum((xt - yt).^2)))/nres );
end

av_rmse = mean(rmses);

end


text:
function varargout=sort(varargin)
%SORT (overloaded)
%
% [t,loc] = sort(x)
%
% The variable t will be the sorted version of x.
%
% SORT is implemented in the nonlinear operator framework using a big-M
% model.

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/SORT CALLED WITH DOUBLE. Report error')

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        if nargin > 1 | min(size(varargin{1}))>1
            error('SDPVAR/SORT only supports simple 1-D sorting'),
        end

        x = varargin{1};
        data.D = binvar(length(x),length(x),'full');
        data.V = sdpvar(length(x),length(x),'full');
        y = [];

        for i = 1:length(x)
            data.i = i;
            data.isthisloc = 0;
            y = [y;yalmip('addextendedvariable',mfilename,x,data)];%i,P,V)];
        end
        loc = [];
        for i = 1:length(x)
            data.i = i;
            data.isthisloc = 1;
            loc = [loc;yalmip('addextendedvariable',mfilename,x,data)];
        end
        varargout{1} = y;
        varargout{2} = loc;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        switch varargin{1}
            case {'milp','graph'}
                % Description using epigraphs
                t = varargin{2};
                X = varargin{3};
                data = varargin{4};
                %i = varargin{4};
                %D = varargin{5};
                %V = varargin{6};

                % Call external to allow subsrefs in classs
                [F,vars] = sort_internal(t,X,data);%i,D,V);

                varargout{1} = F;
                varargout{2} = struct('convexity','milp','monotonicity','milp','definiteness','milp');
                varargout{3} = X;

                % Currently a hack, gneral feature coming soon...
                varargout{2}.models = vars;

            otherwise
                error('SDPVAR/SORT called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/SORT');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function ADMINcreateDictionary(D, HOMEIMAGES, thSize)
% VISUAL OBJECT DICTIONARY
%
% Creates a visual index of all the objects in the database
%
% ADMINcreateDictionary(D, HOMEIMAGES)

homethumbnails = 'objectDictionary';
if nargin<3
    thSize = 64;
    thSize = 48;
end
minNumInstances = 0;

% List of objects
[names, counts, imagendx, objectndx, descriptionndx] = LMobjectnames(D);

% remove unfrequent objects
if minNumInstances>0
    j = find(counts>=minNumInstances);
    names  = names(j);
    counts = counts(j);
    jj = find(ismember(descriptionndx, j));
    imagendx  = imagendx(jj);
    objectndx = objectndx(jj);
    descriptionndx = descriptionndx(jj);
    [foo,descriptionndx] = ismember(descriptionndx, unique(descriptionndx));
end

% sort objects
[foo, ndx] = sort(counts, 'descend');
%[foo, ndx] = sort(names);
Nnames = length(names);
Ninstances = sum(counts);%length(imagendx);

% Create index with global names and with local instances
fid = fopen('dictionary.js', 'w');
fprintf(fid, 'var d=new Array();\n')
fprintf(fid, 'var n=new Array();\n')
fprintf(fid, 'var x=new Array();\n')
fprintf(fid, 'var a=new Array();\n')
fprintf(fid, 'var p=new Array();\n')
fprintf(fid, 'var q=new Array();\n')
fprintf(fid, 'var Nnames=%d;\n', Nnames)
fprintf(fid, 'var Ninstances=%d;\n', Ninstances)
fprintf(fid, 'var nrows=%d;\n', Nnames)
fprintf(fid, 'var ncols=%d;\n', 1)

C = 0;
for i = 1:Nnames
    na = strrep(names{ndx(i)}, '''', '');
    if length(na)>1
        na = [upper(na(1)) lower(na(2:end))];
        fprintf(fid, 'd[%d]=''%s'';\n', i, na);          % strings, object names
        fprintf(fid, 'n[%d]=''%d'';\n', i, counts(ndx(i)));         % number of instances
        fprintf(fid, 'x[%d]=''%d'';\n', i, ndx(i)); % instance index
        fprintf(fid, 'a[%d]=%d;\n', i, C); % instance index
        C = C+counts(ndx(i));
    else
        disp('short name:')
        disp(na)
    end
end

m=0;
for i = 1:Nnames
    k = find(descriptionndx==ndx(i));
    for j= 1:length(k)
        m = m+1;
        o = objectndx(k(j));
        im = imagendx(k(j));
        %page = [D(im).annotation.folder '&image=' D(im).annotation.filename];
        %page = strrep(page, '.jpg', '');
        fprintf(fid, 'p[%d]="%s";\n', m, D(im).annotation.folder);  % strings, object names
        fprintf(fid, 'q[%d]="%s";\n', m, strrep(D(im).annotation.filename, '.jpg', ''));  % strings, object names
    end
end

fclose(fid)

% Loop on D, loop on object
instancesCounts = zeros(Nnames,1);
problems = {''};
for i = 1:length(D)
    disp(i)
    k = find(imagendx==i);

    if ~isempty(k)
        o = objectndx(k);
        d = descriptionndx(k);

        % get thumbnail and extract object and mask
        try
            [annotation, img] = LMread(D, i, HOMEIMAGES);
            for n = 1:length(o)
                %[imgCrop, scaling] = LMobjectnormalizedcrop(img, annotation, o(n), 2, 64, 64);
                imgCrop = LMobjectcrop(img, annotation, o(n), 3);
                %T = imresize(imgCrop, thSize/max(size(imgCrop)), 'bilinear');
                %[foo, T] = LMimpad('', T, [thSize thSize]+6, 255);
                T = imresize(imgCrop, thSize/size(imgCrop,1), 'bilinear');
                [foo, T] = LMimpad('', T, [size(T,1) size(T,2)]+2, 255);

                %thumb = LMobjectThumbnail(D(i).annotation, img, o(n), [64 64]);

                %T = thumb(:,1:end/2-1,:);
                instancesCounts(d(n)) = instancesCounts(d(n))+1;

                %folder = fullfile(homethumbnails, sprintf('obj_%d', d(n)));
                na = strrep(names{d(n)}, '''', '');
                if length(na)>1
                    na = [upper(na(1)) lower(na(2:end))];
                    %folder = fullfile(homethumbnails, strrep(na, ' ', '_'));
                    folder = fullfile(homethumbnails, sprintf('obj_%d', d(n)));

                    filename = fullfile(folder, sprintf('ins_%d.jpg', instancesCounts(d(n))));

                    mkdir(folder)
                    filename
                    imwrite(T, filename, 'jpg', 'quality', 90)

                    %imshow(T)
                    %size(T)
                    %drawnow

                    % Add to description file: folder, SUN name
                    fid = fopen(fullfile(folder, 'description.txt'), 'a');
                    fprintf(fid, '%s/%s \n', D(i).annotation.folder, D(i).annotation.filename)
                    fclose(fid)
                end
            end
        catch
            a = lasterror
            problems = [problems; a.message];
        end
    end
end

% Loop on folders and sort object instances by prototypicality. We can
% estimate prototypicality by training a classifier and sorting the samples
% using the classification score.

author: CSAILVision

repo: LabelMeToolbox


text:
function [params] = kj_calc(params)
% function [params] = kj_calc(params)
% -----------------------------------
% Calculation of the Kagan & Jackson forecasting test.
%
% Input parameters:
%   params.mCatalog           Earthquake catalog
%   params.mPolygon           Polygon (defined by ex_selectgrid)
%   params.vX                 X-vector (defined by ex_selectgrid)
%   params.vY                 Y-vector (defined by ex_selectgrid)
%   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
%   params.bRandom            Perform random simulation (=1) or real calculation (=0)
%   params.nCalculation       Number of random simulations
%   params.bMap               Calculate a map (=1) or a cross-section (=0)
%   params.bNumber            Use constant number (=1) or constant radius (=0)
%   params.nNumberEvents      Number of earthquakes if bNumber == 1
%   params.fRadius            Radius of gridnode if bNumber == 0
%   params.nMinimumNumber     Minimum number of earthquakes per node for determining a b-value
%   params.fForecastPeriod    Forecasting period in years
%   params.bLearningPeriod    Use different learning period than the rest of the catalog
%   params.fLearningPeriod    Learning period in years
%   params.bSignificance      Calculate significance during random simulation
%                             using params.fRealProbability
%   params.fRealProbability   Probability of calculation with real data
%   params.nCalculateMC       Method to calculate magnitude of completeness (see also: help calc_Mc)
%   params.nTestMethod        Method to calculate the Kagan & Jackson test (see also: help kj_poissonian)
%   params.bMinMagMc          Use magnitude of completeness as lower limit of magnitude range for testing (=1)
%                             Use params.fMinMag as lower limit (=0)
%   params.fMinMag            Lower limit of magnitude range for testing
%   params.fMaxMag            Upper limit of magnitude range for testing
%
% Output parameters:
%   Same as input parameters including
%   params.mValueGrid         Matrix of calculated Kagan & Jackson test values
%   params.vRandomMeans       Vector of means of probability differences per simulation run
%   params.vSignificanceLevel Vector of significance levels per simulation run
%   params.fBValueOverall     Calculated overall b-value
%   params.fStdDevOverall     Calculated standard deviation
%   params.fMcOverall         Calculated magnitude of completeness
%
% Danijel Schorlemmer
% March 13, 2002

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Init result matrix
params.mValueGrid = [];
params.vRandomMeans = [];
params.vSignificanceLevel = [];
params.vNormSignificanceLevel = [];

if params.bRandom
  % Create the catalogs for each node with pointers to the overall catalog
  nNumberNodes = length(params.mPolygon(:,1));
  caNodeIndices = cell(nNumberNodes, 1);

  % If cross-section calculate the lenght along cross-section
  if ~params.bMap
    [nRow, nColumn] = size(params.mCatalog);
    xsecx2 = params.mCatalog(:,nColumn);  % length along x-section
    xsecy2 = params.mCatalog(:,7);        % depth of hypocenters
  end

  % loop over all points of the polygon
  disp(['Creating ' num2str(nNumberNodes) ' subcatalogs']);
  for nNode = 1:nNumberNodes

    x = params.mPolygon(nNode, 1);
    y = params.mPolygon(nNode, 2);

    % Calculate distance from center point and sort with distance
    if params.bMap
      vDistances = sqrt(((params.mCatalog(:,1)-x)*cos(pi/180*y)*111).^2 + ((params.mCatalog(:,2)-y)*111).^2);
    else
      vDistances = sqrt(((xsecx2 - x)).^2 + ((xsecy2 + y)).^2);
    end
    if params.bNumber
      % Use first nNumberEvents events
      [vTmp, vIndices] = sort(vDistances);
      caNodeIndices{nNode} = vIndices(1:params.nNumberEvents);
    else
      % Use all events within fRadius
      caNodeIndices{nNode} = find(vDistances <= params.fRadius);
    end
  end % of for nNode
  disp([num2str(nNumberNodes) ' subcatalogs created']);

  % Init values
  vProbDiff = NaN(length(params.mPolygon(:,1)),1);
  vProbK = NaN(length(params.mPolygon(:,1)),1);
  vProbO = NaN(length(params.mPolygon(:,1)),1);

  % overall b-value and Mc
  [v1 params.fBValueOverall params.fStdDevOverall v2] =  bmemag(params.mCatalog);
  params.fMcOverall = calc_Mc(params.mCatalog, params.nCalculateMC);

  % Determine the overall standard deviation of the b-value for the bayesian approach
  if params.nTestMethod == 3
    disp(['Calculating overall standard deviation']);
    params.fStdDevOverall = kj_CalcOverallStdDev(params);
    disp(['Standard deviation calculated']);
  end

  % Define the maximum length of periods
  fMaxLearning = params.fSplitTime - min(params.mCatalog(:,3));
  fMaxForecast = max(params.mCatalog(:,3)) - params.fSplitTime;

  % Do the loop over all calculations
  for nCnt = 1:params.nCalculation
    % Init some variables
    params.mValueGrid = [];
    % Init values
    vProbDiff = NaN(length(params.mPolygon(:,1)),1);
    vProbK = NaN(length(params.mPolygon(:,1)),1);
    vProbO = NaN(length(params.mPolygon(:,1)),1);

    % Permute magnitudes
    mRandomCatalog = params.mCatalog;
    mRandomCatalog(:,6) = params.mCatalog(randperm(length(params.mCatalog)), 6);

    % loop over all points of the polygon
    for nNode = 1:length(params.mPolygon(:,1))

      % Create node catalog
      mNodeCatalog = mRandomCatalog(caNodeIndices{nNode}, :);

      % Determine the local magnitude of completeness
      fMc = calc_Mc(mNodeCatalog, params.nCalculateMC);
      if isnan(fMc)
        fMc = params.fMcOverall;
      elseif isempty(fMc)
        fMc = params.fMcOverall;
      end

      % Create the learning and observed catalogs
      [mLearningCatalog, mObservedCatalog] = ex_SplitCatalog(mNodeCatalog, params.fSplitTime, ...
        params.bLearningPeriod, params.fLearningPeriod, ...
        params.bForecastPeriod, params.fForecastPeriod);

      % Adjust the periods (must not be longer than the catalog contains data)
      if params.bLearningPeriod
        params.fLearning = min(params.fLearningPeriod, fMaxLearning);
      else
        params.fLearning = fMaxLearning;
      end
      if params.bForecastPeriod
        params.fForecast = min(params.fForecastPeriod, fMaxForecast);
      else
        params.fForecast = fMaxForecast;
      end

      % Define magnitude range for testing
      if params.bMinMagMc
        fMinMag = fMc;
      else
        fMinMag = params.fMinMag;
      end

      % Do the Kagan & Jackson test
      [fProbDiff, fProbK, fProbO] = kj_poissonian(mLearningCatalog, params.fLearning, mObservedCatalog, params.fForecast, params.nTestMethod, ...
        params.nMinimumNumber, fMc, params.fBValueOverall, params.fStdDevOverall, fMinMag, params.fMaxMag, 0);
      if ~isnan(fProbDiff)
        if isnan(vProbDiff(nNode))
          vProbDiff(nNode) = 0;
        end
        vProbDiff(nNode) = vProbDiff(nNode) + (fProbDiff/params.nCalculation);
      end
      if ~isnan(fProbK)
        if isnan(vProbK(nNode))
          vProbK(nNode) = 0;
        end
        vProbK(nNode) = vProbK(nNode) + (fProbK/params.nCalculation);
      end
      if ~isnan(fProbO)
        if isnan(vProbO(nNode))
          vProbO(nNode) = 0;
        end
        vProbO(nNode) = vProbO(nNode) + (fProbO/params.nCalculation);
      end
      params.mValueGrid = [params.mValueGrid; vProbDiff(nNode) vProbK(nNode) vProbO(nNode)];
    end  % for nNode
    params.vRandomMeans = [params.vRandomMeans; mean(params.mValueGrid(:, 1), 'omitnan')];
    disp(['Run #' num2str(nCnt) ' of ' num2str(params.nCalculation) ' calculated']);
  end % of for nCalculation
  if params.bSignificance
    nSignificanceLevel = kj_calcsig(params.fRealProbability, params.vRandomMeans, 4);
    disp(['Level of significance: ' num2str(nSignificanceLevel) '% at run #: ' num2str(nCnt) ' Period: ' num2str(params.fForecastPeriod) ' Radius: ' num2str(params.fRadius)]);
    params.vSignificanceLevel = [params.vSignificanceLevel; nSignificanceLevel];
    nSignificanceLevel = kj_CalcNormSig(params.vRandomMeans, params.fRealProbability);
    disp(['Level of significance: ' num2str(nSignificanceLevel) '% at run #: ' num2str(nCnt) ' Period: ' num2str(params.fForecastPeriod) ' Radius: ' num2str(params.fRadius)]);
    params.vNormSignificanceLevel = [params.vNormSignificanceLevel; nSignificanceLevel];
  end
  params.vcsGridNames = cellstr(char('Random: Probability difference (mean)', 'Random: Kagan & Jackson (mean)', 'Random: Our model (mean)'));
else % of if bRandom
  % Determine the overall b-value and magnitude of completeness
  [v1 params.fBValueOverall params.fStdDevOverall v2] =  bmemag(params.mCatalog);
  params.fMcOverall = calc_Mc(params.mCatalog, params.nCalculateMC);

  % Determine the overall standard deviation of the b-value for the bayesian approach (if used)
  if params.nTestMethod == 3
    disp(['Calculating overall standard deviation']);
    params.fStdDevOverall = kj_CalcOverallStdDev(params);
    disp(['Standard deviation calculated']);
  end

  % Define the maximum length of periods
  fMaxLearning = params.fSplitTime - min(params.mCatalog(:,3));
  fMaxForecast = max(params.mCatalog(:,3)) - params.fSplitTime;

  % Loop over all grid nodes
  for nNode = 1:length(params.mPolygon(:,1))
    x = params.mPolygon(nNode, 1);
    y = params.mPolygon(nNode, 2);

    % Calculate distances from center point
    if params.bMap
      vDistances = sqrt(((params.mCatalog(:,1)-x)*cos(pi/180*y)*111).^2 + ((params.mCatalog(:,2)-y)*111).^2);
    else
      [nRow, nColumn] = size(params.mCatalog);
      xsecx2 = params.mCatalog(:,nColumn);                        % Length along cross-section
      xsecy2 = params.mCatalog(:,7);                              % Depth of hypocenters
      vDistances = sqrt(((xsecx2 - x)).^2 + ((xsecy2 + y)).^2);
    end

    % Select the events for calculation
    if params.bNumber
      % Sort with distance
      [vTmp, vIndices] = sort(vDistances);
      mNodeCatalog = params.mCatalog(vIndices(:,1),:);
      % Use first nNumberEvents events
      mNodeCatalog = mNodeCatalog(1:params.nNumberEvents,:);
    else
      % Use all events within fRadius
      vDistances = (vDistances <= params.fRadius);
      mNodeCatalog = params.mCatalog(vDistances,:);
    end

    % Determine the local magnitude of completeness
    fMc = calc_Mc(mNodeCatalog, params.nCalculateMC);
    if isnan(fMc)
      fMc = params.fMcOverall;
    elseif isempty(fMc)
      fMc = params.fMcOverall;
    end

    % Create the learning and observed catalogs
    [mLearningCatalog, mObservedCatalog] = ex_SplitCatalog(mNodeCatalog, params.fSplitTime, ...
      params.bLearningPeriod, params.fLearningPeriod, ...
      params.bForecastPeriod, params.fForecastPeriod);

    % Adjust the periods (must not be longer than the catalog contains data)
    if params.bLearningPeriod
      params.fLearning = min(params.fLearningPeriod, fMaxLearning);
    else
      params.fLearning = fMaxLearning;
    end
    if params.bForecastPeriod
      params.fForecast = min(params.fForecastPeriod, fMaxForecast);
    else
      params.fForecast = fMaxForecast;
    end

    % Define magnitude range for testing
    if params.bMinMagMc
      fMinMag = fMc;
    else
      fMinMag = params.fMinMag;
    end

    % Do the Kagan & Jackson test
    [fProbDiff, fProbK, fProbO, fWeightK, fWeightO, fBValueO] = kj_test(mLearningCatalog, params.fLearning, ...
      mObservedCatalog, params.fForecast, params.nTestMethod, params.nMinimumNumber, fMc, ...
      params.fBValueOverall, params.fStdDevOverall, fMinMag, params.fMaxMag, 0);

    nNumberEventsLearning = length(mLearningCatalog(:,6));
    nNumberEventsObserved = length(mObservedCatalog(:,6));

    % Store the results
    params.mValueGrid = [params.mValueGrid; fProbDiff fProbK fProbO nNumberEventsLearning nNumberEventsObserved fWeightK fWeightO params.fBValueOverall fBValueO];
    params.vcsGridNames = cellstr(char('Probability difference', 'Kagan & Jackson', 'Our model', 'Number events in learning period', ...
      'Number events in forecasting period', 'Weighting of the overall b-value', 'Weighting of the node b-value', ...
      'b-value used for Kagan & Jackson model', 'b-value used for our model'));
  end % for nNode
end

author: CelsoReyes


text:
%% OMAS SYSTEM ENUMS (OMAS_systemType.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef OMAS_systemType < uint8
   enumeration
      windows(0);
      OSX(1);
      linux(2);
   end
end


text:
function V = coeffs2vals(C)
%VAL2COEFFS   Convert tensor of Chebyshev coefficients to values.
%   V = COEFFS2VALS(C) converts a tensor C of bivariate Chebyshev coefficients
%   to a matrix of samples V corresponding to values of
%                   \sum_i \sum_j C(i,j) T_{i-1}(y) T_{j-1}(x)
%   at a tensor Chebyshev grid of size size(C).
%
% See also chebfun3/vals2coeffs.

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

[m, n, p] = size(C);

%% Step 1:
% Mode-1 unfolding of F to get a matrix of size m x n*p:
C1 = chebfun3.unfold(C, 1);

% Apply 1D "vals2coeffs" in the X direction i.e., the 1st dimension of vals1:
C1 = chebtech2.coeffs2vals(C1);

% Tensorize vals1 back to its original m x n x p size:
C1 = chebfun3.fold(C1, [m, n, p], 1, [2 3]);
% This is simply vals1 = reshape(V1, m, n, p);

%% Step 2:
% Mode-2 unfolding of (the tensorized) vals1 to get a matrix of size m x n*p:
C2 = chebfun3.unfold(C1, 2);

% Apply 1D "vals2coeffs" in the Y direction, i.e. to the 1st direction of
% vals2:
C2 = chebtech2.coeffs2vals(C2);

% Tensorize vals2 back to its original m x n x p size:
C2 = chebfun3.fold(C2, [m, n, p], 2, [1 3]);

%% Step 3:
% Mode-3 unfolding of (the tensorized) vals2 to get a matrix of size p x m*n:
C3 = chebfun3.unfold(C2, 3);

% Now, vals2coeffs is applied in the Z direction.
C3 = chebtech2.coeffs2vals(C3);

% Reshape vals3 back to the original m x n x p size:
V = chebfun3.fold(C3, [m, n, p], 3, [1 2]);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function quadrule_test096 ( )

%*****************************************************************************80
%
%% QUADRULE_TEST096 tests HERMITE_GK**_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_TEST096\n' );
  fprintf ( 1, '  HERMITE_GK**_SET sets up a nested rule\n' );
  fprintf ( 1, '  for the Hermite integration problem.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The integration interval is ( -oo, +oo ).\n' );
  fprintf ( 1, '  The weight function is exp ( - x * x )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  HERMITE_INTEGRAL determines the exact value of\n' );
  fprintf ( 1, '  the integal when f(x) = x^m.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here, we just test the highest order rule.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  HERMITE_GK16_SET:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         M         N       Estimate       Exact            Error\n' );
  fprintf ( 1, '\n' );

  n = 35;
  p = 51;

  [ x, w ] = hermite_gk16_set ( n );

  for m = 0 : 2 : min ( p + 2, 20 )

    exact = hermite_integral ( m );

    if ( m == 0 )
      f = ones ( n, 1 );
    else
      f = zeros ( n, 1 );
      f = x.^m;
    end

    estimate = w(1:n)' * f(1:n);

    error = abs ( exact - estimate );

    fprintf ( 1, '  %8d  %8d  %14.6e  %14.6e  %14.6e\n', ...
      m, n, estimate, exact, error );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  HERMITE_GK18_SET:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         M         N       Estimate       Exact            Error\n' );
  fprintf ( 1, '\n' );

  n = 37;
  p = 55;

  [ x, w ] = hermite_gk18_set ( n );

  for m = 0 : 2 : min ( p + 2, 20 )

    exact = hermite_integral ( m );

    if ( m == 0 )
      f = ones ( n, 1 );
    else
      f = zeros ( n, 1 );
      f = x.^m;
    end

    estimate = w(1:n)' * f(1:n);

    error = abs ( exact - estimate );

    fprintf ( 1, '  %8d  %8d  %14.6e  %14.6e  %14.6e\n', ...
      m, n, estimate, exact, error );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  HERMITE_GK22_SET:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         M         N       Estimate       Exact            Error\n' );
  fprintf ( 1, '\n' );

  n = 41;
  p = 63;

  [ x, w ] = hermite_gk22_set ( n );

  for m = 0 : 2 : min ( p + 2, 20 )

    exact = hermite_integral ( m );

    if ( m == 0 )
      f = ones ( n, 1 );
    else
      f = zeros ( n, 1 );
      f = x.^m;
    end

    estimate = w(1:n)' * f(1:n);

    error = abs ( exact - estimate );

    fprintf ( 1, '  %8d  %8d  %14.6e  %14.6e  %14.6e\n', ...
      m, n, estimate, exact, error );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  HERMITE_GK24_SET:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         M         N       Estimate       Exact            Error\n' );
  fprintf ( 1, '\n' );

  n = 43;
  p = 67;

  [ x, w ] = hermite_gk24_set ( n );

  for m = 0 : 2 : min ( p + 2, 20 )

    exact = hermite_integral ( m );

    if ( m == 0 )
      f = ones ( n, 1 );
    else
      f = zeros ( n, 1 );
      f = x.^m;
    end

    estimate = w(1:n)' * f(1:n);

    error = abs ( exact - estimate );

    fprintf ( 1, '  %8d  %8d  %14.6e  %14.6e  %14.6e\n', ...
      m, n, estimate, exact, error );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_bug2741

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_read_cifti ft_write_cifti

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2741'));

cii1 = ft_read_cifti('zstat1.dtseries.nii', 'cortexleft', '100307.L.midthickness.32k_fs_LR.surf.gii', 'cortexright', '100307.R.midthickness.32k_fs_LR.surf.gii');

ft_write_cifti('zstat2', cii1, 'parameter', 'dtseries');

cii2 = ft_read_cifti('zstat2.dtseries.nii'); % surfaces will be read automatically

cii1 = rmfield(cii1, 'hdr');
cii2 = rmfield(cii2, 'hdr');

assert(isequaln(cii1, cii2));