text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef responseClass<handle
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The  ``responseClass`` creates an ``output`` object saved to the MATLAB workspace
    % that contains structures for each instance of a WEC-Sim class (e.g.
    % ``waveClass``, ``bodyClass``, ``constraintClass``, ``ptoClass``,
    % ``cableClass``, ``mooringClass``, etc).
    %
    %.. autoattribute:: objects.responseClass.responseClass
    %
    %.. autoattribute:: objects.responseClass.wave
    %
    %   * ``type`` (`string`) = 'waveType'
    %   * ``time`` (`array`) = [# of time-steps x 1]
    %   * ``elevation`` (`array`) = [# of time-steps x 1]
    %   * ``waveGauge1Elevation`` (`array`) = [# of time-steps x 1] Wave elevation at the location of wave gauge 1
    %   * ``waveGauge2Elevation`` (`array`) = [# of time-steps x 1] Wave elevation at the location of wave gauge 2
    %   * ``waveGauge3Elevation`` (`array`) = [# of time-steps x 1] Wave elevation at the location of wave gauge 3
    %
    %.. autoattribute:: objects.responseClass.bodies
    %
    %   * ``name`` (`string`) = 'bodyName'
    %   * ``time`` (`array`) = [# of time-steps x 1]
    %   * ``position`` (`array`) = [# of time-steps x 6]
    %   * ``velocity`` (`array`) = [# of time-steps x 6]
    %   * ``acceleration`` (`array`) = [# of time-steps x 6]
    %   * ``forceTotal`` (`array`) = [# of time-steps x 6] The sum of all hydrodynamic forces applied to the body
    %   * ``forceExcitation`` (`array`) = [# of time-steps x 6] The sum of the Froude-Krylov excitation force and the mean drift force exerted by the incoming wave on the body
    %   * ``forceRadiationDamping`` (`array`) = [# of time-steps x 6] The negative radiation damping force due to body velocity
    %   * ``forceAddedMass`` (`array`) = [# of time-steps x 6] The negative added mass force due to body acceleration
    %   * ``forceRestoring`` (`array`) = [# of time-steps x 6] The negative sum of the gravity force, buoyant force, the hydrostatic stiffness force, and any moment due to separation between the center of gravity and the center of buoyancy
    %   * ``forceMorisonAndViscous`` (`array`) = [# of time-steps x 6] The negative sum of the Morison element force and the viscous drag force
    %   * ``forceLinearDamping`` (`array`) = [# of time-steps x 6] The negative force due to linear damping and the body velocity
    %
    %   There are 4 additional ``output.bodies`` arrays when using nonlinear hydro and Paraview output:
    %
    %   * ``cellPressures_time`` (`array`) = [# of Paraview time-steps x 1] Nonlinear calculation timeseries
    %   * ``cellPressures_hydrostatic`` (`array`) = [# of Paraview time-steps x # of mesh faces] Hydrostatic pressure on each stl facet
    %   * ``cellPressures_waveLinear`` (`array`) = [# of Paraview time-steps x # of mesh faces] Excitation pressure on each stl facet given zero displacement and the mean free surface
    %   * ``cellPressures_waveNonLinear`` (`array`) = [# of Paraview time-steps x # of mesh faces] Excitation pressure on each stl facet given the instantaneous displacement and instantaneous free surface
    %
    %.. autoattribute:: objects.responseClass.constraints
    %
    %   * ``name`` (`string`) = 'constraintName'
    %   * ``time`` (`array`) = [# of time-steps x 1]
    %   * ``position`` (`array`) = [# of time-steps x 6] The constraint position relative to the initial condition
    %   * ``velocity`` (`array`) = [# of time-steps x 6] The constraint velocity relative to the initial condition
    %   * ``acceleration`` (`array`) = [# of time-steps x 6] The constraint acceleration relative to the initial condition
    %   * ``forceConstraint`` (`array`) = [# of time-steps x 6] The force required to resist motion in the restricted DOFs
    %
    %.. autoattribute:: objects.responseClass.ptos
    %
    %   * ``name`` (`string`) = 'ptoName'
    %   * ``time`` (`array`) = [# of time-steps x 1]
    %   * ``position`` (`array`) = [# of time-steps x 6] The constraint position relative to the initial condition
    %   * ``velocity`` (`array`) = [# of time-steps x 6] The constraint velocity relative to the initial condition
    %   * ``acceleration`` (`array`) = [# of time-steps x 6] The constraint acceleration relative to the initial condition
    %   * ``forceTotal`` (`array`) = [# of time-steps x 6] The sum of the actuation, constraint and internal mechanics forces
    %   * ``forceActuation`` (`array`) = [# of time-steps x 6] The prescribed force input to the PTO to control its motion
    %   * ``forceConstraint`` (`array`) = [# of time-steps x 6] The force required to resist motion in the restricted DOFs
    %   * ``forceInternalMechanics`` (`array`) = [# of time-steps x 6] The net force in the joint DOF due to stiffness and damping
    %   * ``powerInternalMechanics`` (`array`) = [# of time-steps x 6] The net power lost in the joint DOF due to stiffness and damping
    %
    %.. autoattribute:: objects.responseClass.cables
    %
    %   * ``name`` (`string`) = 'cableName'
    %   * ``time`` (`array`) = [# of time-steps x 1]
    %   * ``position`` (`array`) = [# of time-steps x 6]
    %   * ``velocity`` (`array`) = [# of time-steps x 6]
    %   * ``acceleration`` (`array`) = [# of time-steps x 6]
    %   * ``forceTotal`` (`array`) = [# of time-steps x 6] The sum of the actuation and constraint forces
    %   * ``forceActuation`` (`array`) = [# of time-steps x 6] The cable tension
    %   * ``forceConstraint`` (`array`) = [# of time-steps x 6] The force required to resist motion in the restricted DOFs
    %
    %.. autoattribute:: objects.responseClass.mooring
    %
    %   * ``position`` (`array`) = [# of time-steps x 6]
    %   * ``velocity`` (`array`) = [# of time-steps x 6]
    %   * ``forceMooring`` (`array`) = [# of time-steps x 6] The sum of the stiffness, damping and pretension forces applied on the body by the mooring
    %
    %.. autoattribute:: objects.responseClass.moorDyn
    %
    %   * ``Lines`` (`struct`) = [1 x 1] Contains the time and fairlead tensions
    %   * ``Line#`` (`struct`) = [1 x 1] One structure for each mooring line: Line1, Line2, etc. Each line structure contains node positions in x, y, z and segment tensions
    %
    %.. autoattribute:: objects.responseClass.ptosim
    %
    %   * ``time`` (`struct`) = [# of time-steps x 1] Simulation timeseries
    %
    %   There are additional ``output.ptosim`` structs corresponding to the Simulink blocks used:
    %
    %   * ``pistonCF`` (`struct`) = [1 x # of pistons] Structure containing timeseries of compressible fluid piston properties including absolute power, force, position, velocity
    %   * ``pistonNCF`` (`array`) = [1 x # of pistons] Structure containing timeseries of non-compressible fluid piston properties including absolute power, force, top pressure and bottom pressure
    %   * ``checkValve`` (`struct`) = [1 x # of valves] Structure containing timeseries of check valve properies including volume flow rate
    %   * ``valve`` (`struct`) = [1 x # of valves] Structure containing timeseries of valve properties including volume flow rate
    %   * ``accumulator`` (`struct`) = [1 x # of accumulators] Structure containing timeseries of accumulator properties including pressure and volume
    %   * ``hydraulicMotor`` (`struct`) = [1 x # of motors] Structure containing timeseries of hydraulic motor properties including angular velocity and volume flow rate
    %   * ``rotaryGenerator`` (`struct`) = [1 x # of generators] Structure containing timeseries of rotary generator properties including electrical power and generated power
    %   * ``pmLinearGenerator`` (`struct`) = [1 x # of generators] Structure containing timeseries of permanent magnet linear generator properties including absolute power, force, friction force, current, voltage, velocity and electrical power
    %   * ``pmRotaryGenerator`` (`struct`) = [1 x # of generators] Structure containing timeseries of permanent magnet rotary generator properties including absolute power, force, friction force, current, voltage, velocity and electrical power
    %   * ``motionMechanism`` (`struct`) = [1 x # of mechanisms] Structure containing timeseries of motion mechanism properties including PTO torque, angular position and angular velocity
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (SetAccess = 'public', GetAccess = 'public')
        bodies              = struct()     % This property contains a structure for each instance of the ``bodyClass`` (i.e. for each Body block)
        cables              = struct()     % This property contains a structure for each instance of the ``cableClass`` (i.e. for each Cable block)
        constraints         = struct()     % This property contains a structure for each instance of the ``constraintClass`` (i.e. for each Constraint block). Constraint motion is relative from frame F to frame B. Constraint forces act on frame F.
        moorDyn             = struct()     % This property contains a structure for each instance of the ``mooringClass`` using MoorDyn (i.e. for each MoorDyn block)
        mooring             = struct()     % This property contains a structure for each instance of the ``mooringClass`` using the mooring matrix (i.e. for each MooringMatrix block)
        ptos                = struct()     % This property contains a structure for each instance of the ``ptoClass`` (i.e. for each PTO block). PTO motion is relative from frame F to frame B. PTO forces act on frame F.
        ptosim              = struct()     % This property contains a structure for each instance of the ``ptoSimClass`` (i.e. for each PTO-Sim block).
        wave                = struct()     % This property contains a structure for each instance of the ``waveClass``
    end

    methods (Access = 'public')
        function obj = responseClass(bodiesOutput,ptosOutput,constraintsOutput,ptosimOutput,cablesOutput,mooringOutput,waveOutput)
            % This method initializes the ``responseClass``, reads
            % output from each instance of a WEC-Sim class (e.g.
            % ``waveClass``, ``bodyClass``, ``ptoClass``, ``mooringClass``, etc)
            % , and saves the response to an ``output`` object.
            %
            % Returns
            % ------------
            %     output : obj
            %         responseClass object
            %

            % Wave
            obj.wave.type = waveOutput.type;
            obj.wave.time = waveOutput.waveAmpTime(:,1);
            obj.wave.elevation = waveOutput.waveAmpTime(:,2);

            % Bodies
            signals = {'position','velocity','acceleration','forceTotal','forceExcitation','forceRadiationDamping','forceAddedMass','forceRestoring','forceMorisonAndViscous','forceLinearDamping'};
            for ii = 1:length(bodiesOutput)
                obj.bodies(ii).name = bodiesOutput(ii).name;
                obj.bodies(ii).time = bodiesOutput(ii).time;
                obj.bodies(ii).centerGravity   = bodiesOutput(ii).centerGravity;
                for jj = 1:length(signals)
                    obj.bodies(ii).(signals{jj}) = bodiesOutput(ii).signals.values(:, (jj-1)*6+1:(jj-1)*6+6);
                end
                if(bodiesOutput(ii).yaw==1)
                    for t = 1:length(obj.wave.time)
                        % convert kinematic data from global frame to local frame (for use with yaw when yaw displacements may be large).
                        rotMatYaw = eulXYZ2RotMat(0, 0, obj.bodies(ii).position(t,6));
                        rotMatGlobal = eulXYZ2RotMat(obj.bodies(ii).position(t,4), obj.bodies(ii).position(t,5), obj.bodies(ii).position(t,6));
                        rotMatLocal = rotMatYaw.' * rotMatGlobal;
                        [phiLoc, thetaLoc, psiLoc] = rotMatXYZ2Eul(rotMatLocal); % get orientation in local frame
                        % position in local frame
                        obj.bodies(ii).positionLocal(t,1:3) = rotMatYaw.'*obj.bodies(ii).position(t,1:3).'; % rotate linear position vector from global to local frame
                        obj.bodies(ii).positionLocal(t,4:6) = [phiLoc, thetaLoc, psiLoc];
                        % velocity in local frame
                        obj.bodies(ii).velocityLocal(t,1:3) = rotMatYaw.'*obj.bodies(ii).velocity(t,1:3).';
                        obj.bodies(ii).velocityLocal(t,4:6) = rotMatYaw.'*obj.bodies(ii).velocity(t,4:6).';
                        % acceleration in local frame
                        obj.bodies(ii).accelerationLocal(t,1:3) = rotMatYaw.'*obj.bodies(ii).acceleration(t,1:3).';
                        obj.bodies(ii).accelerationLocal(t,4:6) = rotMatYaw.'*obj.bodies(ii).acceleration(t,4:6).';
                    end
                end
                if ~isempty(bodiesOutput(ii).hspressure)
                    obj.bodies(ii).cellPressures_time = bodiesOutput(ii).hspressure.time;
                    obj.bodies(ii).cellPressures_hydrostatic   = bodiesOutput(ii).hspressure.signals.values;
                    obj.bodies(ii).cellPressures_waveLinear    = bodiesOutput(ii).wpressurel.signals.values;
                    obj.bodies(ii).cellPressures_waveNonLinear = bodiesOutput(ii).wpressurenl.signals.values;
                else
                    obj.bodies(ii).cellPressures_time = [];
                    obj.bodies(ii).cellPressures_hydrostatic   = [];
                    obj.bodies(ii).cellPressures_waveLinear    = [];
                    obj.bodies(ii).cellPressures_waveNonLinear = [];
                end
            end
            % PTOs
            if isstruct(ptosOutput)
                signals = {'position','velocity','acceleration','forceTotal','forceActuation','forceConstraint','forceInternalMechanics','powerInternalMechanics'};
                for ii = 1:length(ptosOutput)
                    obj.ptos(ii).name = ptosOutput(ii).name;
                    obj.ptos(ii).time = ptosOutput(ii).time;
                    for jj = 1:length(signals)
                        obj.ptos(ii).(signals{jj}) =  ptosOutput(ii).signals.values(:,(jj-1)*6+1:(jj-1)*6+6);
                    end
                end
            end
            % Constraints
            if isstruct(constraintsOutput)
                signals = {'position','velocity','acceleration','forceConstraint'};
                for ii = 1:length(constraintsOutput)
                    obj.constraints(ii).name = constraintsOutput(ii).name;
                    obj.constraints(ii).time = constraintsOutput(ii).time;
                    for jj = 1:length(signals)
                        obj.constraints(ii).(signals{jj}) = constraintsOutput(ii).signals.values(:, (jj-1)*6+1:(jj-1)*6+6);
                    end
                end
            end
            % Mooring
            if isstruct(mooringOutput)
                signals = {'position','velocity','forceMooring'};
                for ii = 1:length(mooringOutput)
                    obj.mooring(ii).name = mooringOutput(ii).name;
                    obj.mooring(ii).time = mooringOutput(ii).time;
                    for jj = 1:length(signals)
                        obj.mooring(ii).(signals{jj}) = mooringOutput(ii).signals.values(:, (jj-1)*6+1:(jj-1)*6+6);
                    end
                end
            end
            % Cables
            if isstruct(cablesOutput)
                signals = {'position','velocity','acceleration','forceTotal','forceActuation',...
                    'forceConstraint'};
                for ii = 1:length(cablesOutput)
                    obj.cables(ii).name = cablesOutput(ii).name;
                    obj.cables(ii).time = cablesOutput(ii).time;
                    for jj = 1:length(signals)
                        obj.cables(ii).(signals{jj}) = cablesOutput(ii).signals.values(:,(jj-1)*6+1:(jj-1)*6+6);
                    end
                end
            end
            % PTO-Sim
            if isstruct(ptosimOutput)
                %names = {'HydPistonCompressible','GasHydAccumulator','RectifyingCheckValve','HydraulicMotorV2','ElectricMachineEC'};
                hydPistonCompressibleSignals = {'pressureA','forcePTO','pressureB'};
                gasHydAccumulatorSignals = {'pressure','flowRate'};
                rectifyingCheckValveSignals = {'flowRateA','flowRateB','flowRateC','flowRateD'};
                hydraulicMotorSignals = {'shaftSpeed','torque','deltaP','flowRate'};
                electricMachineECSignals = {'torqueEM','shaftSpeed','current','voltage'};
                linearCrankSignals = {'ptoTorque','angPosition','angVelocity'};
                adjustableRodSignals = {'ptoTorque','angPosition','angVelocity'};
                checkValveSignals = {'flowCheckValve','deltaPCheckValve'};
                linearGeneratorSignals = {'absPower','force','fricForce','Ia','Ib','Ic','Va','Vb','Vc','elecPower','vel'};
                rotaryGeneratorSignals = {'absPower','Torque','fricTorque','Ia','Ib','Ic','Va','Vb','Vc','elecPower','vel'};

                for ii = 1:length(ptosimOutput)
                    obj.ptosim(ii).name = ptosimOutput(ii).name;
                    obj.ptosim(ii).time = ptosimOutput(ii).time;
                    obj.ptosim(ii).type = ptosimOutput(ii).type;
                    if ptosimOutput(ii).type == 1
                        signals = electricMachineECSignals;
                    elseif ptosimOutput(ii).type == 2
                        signals = hydPistonCompressibleSignals;
                    elseif ptosimOutput(ii).type == 3
                        signals = gasHydAccumulatorSignals;
                    elseif ptosimOutput(ii).type == 4
                        signals = rectifyingCheckValveSignals;
                    elseif ptosimOutput(ii).type == 5
                        signals = hydraulicMotorSignals;
                    elseif ptosimOutput(ii).type == 6
                        signals = linearCrankSignals;
                    elseif ptosimOutput(ii).type == 7
                        signals = adjustableRodSignals;
                    elseif ptosimOutput(ii).type == 8
                        signals = checkValveSignals;
                    elseif ptosimOutput(ii).type == 9
                        signals = linearGeneratorSignals;
                    elseif ptosimOutput(ii).type == 10
                        signals = rotaryGeneratorSignals;
                    end
                    for jj = 1:length(signals)
                        obj.ptosim(ii).(signals{jj}) = ptosimOutput(ii).signals.values(:,jj);
                    end
                end
            end
        end

        function obj = loadMoorDyn(obj,linesNum)
            % This method reads MoorDyn outputs for each instance of the
            % ``mooringClass``
            %
            % Parameters
            % ------------
            %     linesNum : integer
            %         the number of MoorDyn lines
            %

            arguments
                obj
                linesNum (1,1) double {mustBeInteger}
            end

            % load Lines.out
            filename = './Mooring/Lines.out';
            fid = fopen(filename, 'r');
            header = strsplit(fgetl(fid));
            data = dlmread(filename,'',1,0);
            tmp = size(data);
            ncol = tmp(2);clear tmp
            for icol=1:ncol
               eval(['obj.moorDyn.Lines.' header{icol} ' = data(:,' num2str(icol) ');']);
            end
            fclose(fid);
            % load Line#.out
            for iline=1:linesNum
                eval(['obj.moorDyn.Line' num2str(iline) '=struct();']);
                filename = ['./Mooring/Line' num2str(iline) '.out'];
                try
                    fid = fopen(filename);
                    header = strsplit(strtrim(fgetl(fid)));
                    data = dlmread(filename,'',1,0);
                    tmp = size(data);
                    ncol = tmp(2);clear tmp
                    for icol=1:ncol
                        eval(['obj.moorDyn.Line' num2str(iline) '.' header{icol} ' = data(:,' num2str(icol) ');']);
                    end
                    fclose(fid);
                catch
                    fprintf('\n No moorDyn *.out file saved for Line%u\n',iline);
                end
            end
        end

        function plotResponse(obj,bodyNum,comp)
            % This method plots the response of a body for a given DOF.
            %
            % Parameters
            % ------------
            %     bodyNum : integer
            %         the body number to plot
            %
            %     comp : integer
            %         the response component (i.e. dof) to be plotted (e.g. 1-6)
            %

            arguments
                obj
                bodyNum (1,1) double {mustBeInteger}
                comp (1,1) double {mustBeMember(comp,1:6)}
            end

            DOF = {'Surge','Sway','Heave','Roll','Pitch','Yaw'};
            t=obj.bodies(bodyNum).time;
            if comp < 4
                pos=obj.bodies(bodyNum).position(:,comp)-obj.bodies(bodyNum).centerGravity(comp);
            else
                pos=obj.bodies(bodyNum).position(:,comp);
            end
            vel=obj.bodies(bodyNum).velocity(:,comp);
            acc=obj.bodies(bodyNum).acceleration(:,comp);
            figure()
            plot(t,pos,'k-',...
                t,vel,'b-',...
                t,acc,'r-')
            legend({'position','velocity','acceleration'})
            xlabel('Time (s)')
            ylabel('Response in (m) or (radians)')
            title(['body' num2str(bodyNum) ' (' obj.bodies(bodyNum).name ') ' DOF{comp} ' Response'])
            clear t pos vel acc i
        end

        function plotForces(obj,bodyNum,comp)
            % This method plots the forces on a body for a given DOF.
            %
            % Parameters
            % ------------
            %     bodyNum : integer
            %         the body number to plot
            %
            %     comp : integer
            %         the force component (i.e. dof) to be plotted (e.g. 1-6)
            %

            arguments
                obj
                bodyNum (1,1) double {mustBeInteger}
                comp (1,1) double {mustBeMember(comp,1:6)}
            end

            DOF = {'Surge','Sway','Heave','Roll','Pitch','Yaw'};
            t=obj.bodies(bodyNum).time;
            FT=obj.bodies(bodyNum).forceTotal(:,comp);
            FE=obj.bodies(bodyNum).forceExcitation(:,comp);
            FRD=-1*obj.bodies(bodyNum).forceRadiationDamping(:,comp);
            FAM=-1*obj.bodies(bodyNum).forceAddedMass(:,comp);
            FR=-1*obj.bodies(bodyNum).forceRestoring(:,comp);
            FMV=-1*obj.bodies(bodyNum).forceMorisonAndViscous(:,comp);
            FLD=-1*obj.bodies(bodyNum).forceLinearDamping(:,comp);
            figure();
            plot(t,FT,...
                t,FE,...
                t,FRD,...
                t,FAM,...
                t,FR,...
                t,FMV,...
                t,FLD)
            legend('forceTotal','forceExcitation','forceRadiationDamping','forceAddedMass','forceRestoring','forceViscous','forceLinearDamping')
            xlabel('Time (s)')
            ylabel('Force (N) or Torque (N*m)')
            title(['body' num2str(bodyNum) ' (' obj.bodies(bodyNum).name ') ' DOF{comp} ' Forces'])
            clear t FT FE FRD FR FV FM i
        end

        function saveViz(obj,simu,body,waves,options)
            % This method plots and saves the wave elevation and body
            % geometry over time to visualize the waves and response
            %
            % Parameters
            % ------------
            %     simu :
            %         simulationClass object
            %
            %     body :
            %         bodyClass object
            %
            %     waves :
            %         waveClass object
            %
            %     options
            %         axisLimits : 1x6 float matrix
            %             x, y, and z-bounds of figure axes in the format:
            %             [min x, max x, min y, max y, min z, max z]
            %             Default = [-simu.domainSize/2 simu.domainSize/2
            %             -simu.domainSize/2 simu.domainSize/2 -waves.waterDepth maximumHeight]
            %         timesPerFrame : float
            %             number of simulation time steps per video frame
            %             (higher number decreases computation time)
            %             Default = 1
            %         startEndTime : 1x2 float matrix
            %             Array defining the start and end times of the
            %             visualization
            %             Default = [min(t) max(t)]
            %         saveSetting : integer
            %             0 = video, 1 = gif. Default = 0
            %

            arguments
                obj
                simu
                body
                waves
                options.axisLimits (1,6) double {mustBeReal, mustBeNonNan, mustBeFinite} = [-simu.domainSize/2 simu.domainSize/2 -simu.domainSize/2 simu.domainSize/2 -waves.waterDepth -999];
                options.timesPerFrame (1,1) double {mustBeInteger, mustBePositive} = 1;
                options.startEndTime (1,2) double {mustBeReal, mustBeNonnegative, mustBeNonNan} = [0 0];
                options.saveSetting (1,1) {mustBeMember(options.saveSetting,0:1)} = 0;
            end
            % Set time vector
            t = obj.wave.time(1:options.timesPerFrame*round(simu.dtOut/simu.dt,0):end,1);
            if isequal(options.startEndTime, [0 0])
                options.startEndTime = [min(t) max(t)];
            end
            % Create grid using provided x and y coordinates
            x = linspace(options.axisLimits(1),options.axisLimits(2),200);
            y = linspace(options.axisLimits(3),options.axisLimits(4),200);
            [X,Y] = meshgrid(x,y);
            % Initialize maximum body height
            maxHeight = 0;
            % Read in data for each body
            for ibod=1:length(obj.bodies)
                % Read and assign geometry data
                if isempty(body(ibod).geometry.vertex)
                    % Import body geometry for linear hydro cases
                    body(ibod).importBodyGeometry(simu.domainSize);
                end
                bodyMesh(ibod).Points = body(ibod).geometry.vertex;
                bodyMesh(ibod).Conns = body(ibod).geometry.face;
                % Read changes and assign angles and position changes over time
                bodyMesh(ibod).deltaPos = [obj.bodies(ibod).position(1:options.timesPerFrame:end,1)-obj.bodies(ibod).position(1,1),...
                obj.bodies(ibod).position(1:options.timesPerFrame:end,2)-obj.bodies(ibod).position(1,2),...
                obj.bodies(ibod).position(1:options.timesPerFrame:end,3)-obj.bodies(ibod).position(1,3)];
                % Save maximum body height
                if max(bodyMesh(ibod).Points(:,3))+max(bodyMesh(ibod).deltaPos(:,3)) > maxHeight
                    maxHeight = max(bodyMesh(ibod).Points(:,3))+max(bodyMesh(ibod).deltaPos(:,3));
                end
            end
            if options.axisLimits(6) == -999
                options.axisLimits(6) = maxHeight;
            end
            if options.saveSetting == 0
                % Create video file and open it for writing
                video = VideoWriter('waveViz.avi');
                video.FrameRate = 1/(simu.dtOut*options.timesPerFrame);
                open(video);
            elseif options.saveSetting == 1
                % Create the gif file
                gifFilename = 'waveViz.gif';
            end
            % Initialize figure and image counter
            figure();
            imageCount = 0;
            for i=1:length(t)
                if t(i) >= options.startEndTime(1) && t(i) <= options.startEndTime(2)
                    imageCount = imageCount + 1;
                    for ibod = 1:length(obj.bodies)
                        % Apply rotation to each point
                        rotMat = eulXYZ2RotMat(obj.bodies(ibod).position(1+options.timesPerFrame*(i-1),4), ...
                            obj.bodies(ibod).position(1+options.timesPerFrame*(i-1),5), ...
                            obj.bodies(ibod).position(1+options.timesPerFrame*(i-1),6));
                        for ipts=1:length(bodyMesh(ibod).Points(:,1))
                            bodyMesh(ibod).rotation(ipts,:) = (rotMat*bodyMesh(ibod).Points(ipts,:).').';
                        end

                        % Calculate full position changes due to rotation,
                        % translation, and center of gravity
                        bodyMesh(ibod).pointsNew = bodyMesh(ibod).rotation + bodyMesh(ibod).deltaPos(i,:) + body(ibod).centerGravity.';

                        % Create and plot final triangulation of geometry with applied changes
                        bodFinal = triangulation(bodyMesh(ibod).Conns,bodyMesh(ibod).pointsNew);
                        trisurf(bodFinal,'FaceColor',[1 1 0],'EdgeColor','k','EdgeAlpha',.2)
                        hold on
                    end
                    % Create and wave elevation grid
                    Z = waveElevationGrid(waves, t(i), X, Y, simu.dtOut, simu.gravity);
                    surf(X,Y,Z,'FaceAlpha',.85,'EdgeColor','none')
                    hold on
                    % Display seafloor
                    seaFloor = -waves.waterDepth*ones(size(X, 1));
                    surf(X,Y,seaFloor,'FaceColor',[.4 .4 0],'EdgeColor','none');
                    hold on
                    % Time visual
                    nDecimals = max(0,ceil(-log10(simu.dtOut*options.timesPerFrame)));
                    nLeading = ceil(log10(max(t)));
                    tAnnot = sprintf(['time = %' num2str(nDecimals+nLeading+1) '.' num2str(nDecimals) 'f s'],t(i));
                    % Settings and labels
                    caxis([min(waves.waveAmpTime(:,2)) max(waves.waveAmpTime(:,2))])
                    colormap winter
                    c = colorbar;
                    ylabel(c, 'Wave Elevation (m)')
                    title({'Wave Elevation and Geometry Visualization',tAnnot})
                    xlabel('x(m)')
                    ylabel('y(m)')
                    zlabel('z(m)')
                    daspect([1 1 1])
                    axis(options.axisLimits)
                    % Create figure while iterating through time loop
                    drawnow;
                    % Capture figure for saving
                    frame = getframe(gcf);
                    if options.saveSetting == 0
                        % Save to video
                        writeVideo(video,frame);
                    elseif options.saveSetting == 1
                        % Save to gif
                        im = frame2im(frame);
                        [imind,cm] = rgb2ind(im,256);
                        if imageCount == 1
                            imwrite(imind,cm,gifFilename,'gif', 'Loopcount',inf);
                        else
                            imwrite(imind,cm,gifFilename,'gif','WriteMode','append','DelayTime',simu.dtOut);
                        end
                    end
                    hold off
                end
            end
            % Close video file
            if options.saveSetting == 0
                close(video);
            end
        end

        function writeText(obj)
            % This method writes WEC-Sim outputs to a (ASCII) text file.
            % This method is executed by specifying ``simu.outputtxt=1``
            % in the ``wecSimInputFile.m``.
            filename = ['output/wave.txt'];
            fid = fopen(filename,'w+');
            header = {'time','elevation'};
            for ii=1:length(header)
                tmp(ii) = length(header{ii});
            end
            numChar = max(tmp)+2; clear tmp;
            header_fmt = ['%' num2str(numChar) 's '];
            ncols = length(header);
            tmp = size(obj.wave.time);
            nrows = tmp(1); clear tmp;
            data_fmt = [repmat('%10.5f ',1,ncols) '\n'];
            data = zeros(nrows,ncols);
            data(:,1) = obj.wave.time;
            data(:,2) = obj.wave.elevation;
            if ncols > 2
                for ii = 3:length(header)
                    eval(['data(:,' num2str(ii) ') = obj.wave.' header{ii} ';']);
                end
            end
            for ii = 1:length(header)
                fprintf(fid,header_fmt,header{ii});
            end
            fprintf(fid,'\n');
            fprintf(fid,data_fmt,data');
            fclose(fid);
            % bodies
            signals = {'position','velocity','acceleration','forceTotal','forceExcitation','forceRadiationDamping','forceAddedMass','forceRestoring','forceMorisonAndViscous','forceLinearDamping'};
            header = {'time', ...
                      'position_1'               ,'position_2'               ,'position_3'               ,'position_4'               ,'position_5'               ,'position_6'               , ...
                      'velocity_1'               ,'velocity_2'               ,'velocity_3'               ,'velocity_4'               ,'velocity_5'               ,'velocity_6'               , ...
                      'acceleration_1'           ,'acceleration_2'           ,'acceleration_3'           ,'acceleration_4'           ,'acceleration_5'           ,'acceleration_6'           , ...
                      'forceTotal_1'             ,'forceTotal_2'             ,'forceTotal_3'             ,'forceTotal_4'             ,'forceTotal_5'             ,'forceTotal_6'             , ...
                      'forceExcitation_1'        ,'forceExcitation_2'        ,'forceExcitation_3'        ,'forceExcitation_4'        ,'forceExcitation_5'        ,'forceExcitation_6'        , ...
                      'forceRadiationDamping_1'  ,'forceRadiationDamping_2'  ,'forceRadiationDamping_3'  ,'forceRadiationDamping_4'  ,'forceRadiationDamping_5'  ,'forceRadiationDamping_6'  , ...
                      'forceAddedMass_1'         ,'forceAddedMass_2'         ,'forceAddedMass_3'         ,'forceAddedMass_4'         ,'forceAddedMass_5'         ,'forceAddedMass_6'         , ...
                      'forceRestoring_1'         ,'forceRestoring_2'         ,'forceRestoring_3'         ,'forceRestoring_4'         ,'forceRestoring_5'         ,'forceRestoring_6'         , ...
                      'forceMorisonAndViscous_1' ,'forceMorisonAndViscous_2' ,'forceMorisonAndViscous_3' ,'forceMorisonAndViscous_4' ,'forceMorisonAndViscous_5' ,'forceMorisonAndViscous_6' , ...
                      'forceLinearDamping_1'     ,'forceLinearDamping_2'     ,'forceLinearDamping_3'     ,'forceLinearDamping_4'     ,'forceLinearDamping_5'     ,'forceLinearDamping_6'     };
            for ii=1:length(signals)
                tmp(ii) = length(signals{ii});
            end
            numChar = max(tmp)+2; clear tmp;
            ncols = 1 + length(signals)*6;
            tmp = size(obj.bodies(1).time);
            nrows = tmp(1); clear tmp;
            header_fmt = ['%' num2str(numChar) 's '];
            data_fmt = [repmat('%10.5f ',1,ncols) '\n'];
            for ibod = 1:length(obj.bodies)
                filename = ['output/body' num2str(ibod) '_' obj.bodies(ibod).name '.txt'];
                fid = fopen(filename,'w+');
                data = zeros(nrows,ncols);
                data(:,1) = obj.bodies(ibod).time;
                fprintf(fid,header_fmt,header{1});
                for isignal=1:length(signals)
                    for idof = 1:6
                        fprintf(fid,header_fmt,header{1 + (isignal-1)*6+idof});
                    end
                    data(:, 1+(isignal-1)*6+1:1+(isignal-1)*6+6) = obj.bodies(ibod).(signals{isignal});
                end
                fprintf(fid,'\n');
                fprintf(fid,data_fmt,data');
                fclose(fid);
                if isfield(obj.bodies(ibod),'cellPressures_hydrostatic')
                    filename = ['output/body' num2str(ibod) '_' obj.bodies(ibod).name '_cellPressure_hydrostatic.txt'];
                    fid = fopen(filename,'w+');
                    header_2 = {'time'};
                    tmp = size(obj.bodies(ibod).cellPressures_hydrostatic);
                    nrows2 = tmp(1);
                    ncols2 = tmp(2)+1;
                    for icell=1:ncols2-1
                        header_2{icell+1} = ['cell_' num2str(icell)];
                    end
                    tmp = length(header_2{end});
                    numChar = max(tmp)+2; clear tmp;
                    header_fmt_2 = ['%' num2str(numChar) 's '];
                    data_fmt_2 = [repmat('%10.5f ',1,ncols2) '\n'];
                    data = zeros(nrows2,ncols2);
                    data(:,1) = obj.bodies(ibod).cellPressures_time;
                    data(:,2:end) = obj.bodies(ibod).cellPressures_hydrostatic;
                    for jj = 1:length(header_2)
                        fprintf(fid,header_fmt_2,header_2{jj});
                    end
                    fprintf(fid,'\n');
                    fprintf(fid,data_fmt_2,data');
                    fclose(fid);
                    % wave linear
                    filename = ['output/body' num2str(ibod) '_' obj.bodies(ibod).name '_cellPressure_waveLinear.txt'];
                    fid = fopen(filename,'w+');
                    data = zeros(nrows2,ncols2);
                    data(:,1) = obj.bodies(ibod).cellPressures_time;
                    data(:,2:end) = obj.bodies(ibod).cellPressures_waveLinear;
                    for jj = 1:length(header_2)
                        fprintf(fid,header_fmt_2,header_2{jj});
                    end
                    fprintf(fid,'\n');
                    fprintf(fid,data_fmt_2,data');
                    fclose(fid);
                    % wave nonlinear
                    filename = ['output/body' num2str(ibod) '_' obj.bodies(ibod).name '_cellPressure_waveNonLinear.txt'];
                    fid = fopen(filename,'w+');
                    data = zeros(nrows2,ncols2);
                    data(:,1) = obj.bodies(ibod).cellPressures_time;
                    data(:,2:end) = obj.bodies(ibod).cellPressures_waveNonLinear;
                    for jj = 1:length(header_2)
                        fprintf(fid,header_fmt_2,header_2{jj});
                    end
                    fprintf(fid,'\n');
                    fprintf(fid,data_fmt_2,data');
                    fclose(fid);
                end
            end
            % ptos
            if isfield(obj.ptos,'name')
                signals = {'position','velocity','acceleration','forceTotal','forceActuation','forceConstraint','forceInternalMechanics','powerInternalMechanics'};
                header = {'time', ...
                          'position_1'              ,'position_2'              ,'position_3'              ,'position_4'              ,'position_5'              ,'position_6'              , ...
                          'velocity_1'              ,'velocity_2'              ,'velocity_3'              ,'velocity_4'              ,'velocity_5'              ,'velocity_6'              , ...
                          'acceleration_1'          ,'acceleration_2'          ,'acceleration_3'          ,'acceleration_4'          ,'acceleration_5'          ,'acceleration_6'          , ...
                          'forceTotal_1'            ,'forceTotal_2'            ,'forceTotal_3'            ,'forceTotal_4'            ,'forceTotal_5'            ,'forceTotal_6'            , ...
                          'forceActuation_1'        ,'forceActuation_2'        ,'forceActuation_3'        ,'forceActuation_4'        ,'forceActuation_5'        ,'forceActuation_6'        , ...
                          'forceConstraint_1'       ,'forceConstraint_2'       ,'forceConstraint_3'       ,'forceConstraint_4'       ,'forceConstraint_5'       ,'forceConstraint_6'       , ...
                          'forceInternalMechanics_1','forceInternalMechanics_2','forceInternalMechanics_3','forceInternalMechanics_4','forceInternalMechanics_5','forceInternalMechanics_6', ...
                          'powerInternalMechanics_1','powerInternalMechanics_2','powerInternalMechanics_3','powerInternalMechanics_4','powerInternalMechanics_5','powerInternalMechanics_6', };
                for ii=1:length(signals)
                    tmp(ii) = length(signals{ii});
                end
                numChar = max(tmp)+2; clear tmp;
                ncols = 1 + length(signals)*6;
                tmp = size(obj.ptos(1).time);
                nrows = tmp(1); clear tmp;
                header_fmt = ['%' num2str(numChar) 's '];
                data_fmt = [repmat('%10.5f ',1,ncols) '\n'];
                for ipto = 1:length(obj.ptos)
                    filename = ['output/pto' num2str(ipto) '_' obj.ptos(ipto).name '.txt'];
                    fid = fopen(filename,'w+');
                    data = zeros(nrows,ncols);
                    data(:,1) = obj.ptos(ipto).time;
                    fprintf(fid,header_fmt,header{1});
                    for isignal=1:length(signals)
                        for idof = 1:6
                            fprintf(fid,header_fmt,header{1 + (isignal-1)*6+idof});
                        end
                        data(:, 1+(isignal-1)*6+1:1+(isignal-1)*6+6) = obj.ptos(ipto).(signals{isignal});
                    end
                    fprintf(fid,'\n');
                    fprintf(fid,data_fmt,data');
                    fclose(fid);
                end
            end
            % constraints
            if isfield(obj.constraints,'name')
                signals = {'position','velocity','acceleration','forceConstraint'};
                header = {'time', ...
                          'position_1'       ,'position_2'       ,'position_3'       ,'position_4'       ,'position_5'       ,'position_6'       , ...
                          'velocity_1'       ,'velocity_2'       ,'velocity_3'       ,'velocity_4'       ,'velocity_5'       ,'velocity_6'       , ...
                          'acceleration_1'   ,'acceleration_2'   ,'acceleration_3'   ,'acceleration_4'   ,'acceleration_5'   ,'acceleration_6'   , ...
                          'forceConstraint_1','forceConstraint_2','forceConstraint_3','forceConstraint_4','forceConstraint_5','forceConstraint_6'     };
                for ii=1:length(signals)
                    tmp(ii) = length(signals{ii});
                end
                numChar = max(tmp)+2; clear tmp;
                ncols = 1 + length(signals)*6;
                tmp = size(obj.constraints(1).time);
                nrows = tmp(1); clear tmp;
                header_fmt = ['%' num2str(numChar) 's '];
                data_fmt = [repmat('%10.5f ',1,ncols) '\n'];
                for icon = 1:length(obj.constraints)
                    filename = ['output/constraint' num2str(icon) '_' obj.constraints(icon).name '.txt'];
                    fid = fopen(filename,'w+');
                    data = zeros(nrows,ncols);
                    data(:,1) = obj.constraints(icon).time;
                    fprintf(fid,header_fmt,header{1});
                    for isignal=1:length(signals)
                        for idof = 1:6
                            fprintf(fid,header_fmt,header{1 + (isignal-1)*6+idof});
                        end
                        data(:, 1+(isignal-1)*6+1:1+(isignal-1)*6+6) = obj.constraints(icon).(signals{isignal});
                    end
                    fprintf(fid,'\n');
                    fprintf(fid,data_fmt,data');
                    fclose(fid);
                end
            end
            %ptoSim
            if isfield(obj.ptosim,'time')
                f1 = fields(obj.ptosim);
                count = 1;
                header = {'time'};
                data = obj.ptosim.time;
                for ifld1=1:(length(f1)-1)
                    f2 = fields(obj.ptosim.(f1{ifld1}));
                    for iins = 1:length(obj.ptosim.(f1{ifld1}))
                        for ifld2 = 1:length(f2)
                            count = count+1;
                            header{count} = [f1{ifld1} num2str(iins) '_' f2{ifld2}];
                            data(:,count) = obj.ptosim.(f1{ifld1}).(f2{ifld2});
                        end
                    end
                end
                for ii=1:length(header)
                    tmp(ii) = length(header{ii});
                end
                numChar = max(tmp)+2; clear tmp;
                header_fmt = ['%' num2str(numChar) 's '];
                data_fmt = [repmat('%10.5f ',1,length(header)) '\n'];
                filename = ['output/ptosim.txt'];
                fid = fopen(filename,'w+');
                for ii=1:length(header)
                    fprintf(fid,header_fmt,header{ii});
                end
                fprintf(fid,'\n');
                fprintf(fid,data_fmt,data');
                fclose(fid);
            end
            % mooring
            if isfield(obj.mooring,'name')
                signals = {'position','velocity','forceMooring'};
                header = {'time', ...
                          'position_1'       ,'position_2'       ,'position_3'       ,'position_4'       ,'position_5'       ,'position_6'       , ...
                          'velocity_1'       ,'velocity_2'       ,'velocity_3'       ,'velocity_4'       ,'velocity_5'       ,'velocity_6'       , ...
                          'forceMooring_1','forceMooring_2','forceMooring_3','forceMooring_4','forceMooring_5','forceMooring_6'     };
                for ii=1:length(signals)
                    tmp(ii) = length(signals{ii});
                end
                numChar = max(tmp)+2; clear tmp;
                ncols = 1 + length(signals)*6;
                tmp = size(obj.mooring(1).time);
                nrows = tmp(1); clear tmp;
                header_fmt = ['%' num2str(numChar) 's '];
                data_fmt = [repmat('%10.5f ',1,ncols) '\n'];
                for imoor = 1:length(obj.mooring)
                    filename = ['output/mooring' num2str(imoor) '_' obj.mooring(imoor).name '.txt'];
                    fid = fopen(filename,'w+');
                    data = zeros(nrows,ncols);
                    data(:,1) = obj.mooring(imoor).time;
                    fprintf(fid,header_fmt,header{1});
                    for isignal=1:length(signals)
                        for idof = 1:6
                            fprintf(fid,header_fmt,header{1 + (isignal-1)*6+idof});
                        end
                        data(:, 1+(isignal-1)*6+1:1+(isignal-1)*6+6) = obj.mooring(imoor).(signals{isignal});
                    end
                    fprintf(fid,'\n');
                    fprintf(fid,data_fmt,data');
                    fclose(fid);
                end
            end
        end
    end
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
% (non-blind) image deblurring

% @inproceedings{zhang2017learning,
%   title={Learning Deep CNN Denoiser Prior for Image Restoration},
%   author={Zhang, Kai and Zuo, Wangmeng and Gu, Shuhang and Zhang, Lei},
%   booktitle={IEEE Conference on Computer Vision and Pattern Recognition},
%   year={2017}
% }

% If you have any question, please feel free to contact with me.
% Kai Zhang (e-mail: cskaizhang@gmail.com)

clear; clc;

addpath('utilities');
imageSets    = {'Set3G','Set3C'}; %%% testing dataset
setTest      = imageSets([1]); %%% select the dataset

showResult   = 1;
pauseTime    = 1;
useGPU       = 1;

folderTest   = 'testsets';
folderResult = 'results';
folderModel  = 'models';
taskTestCur  = 'Deblur';
if ~exist(folderResult,'file')
    mkdir(folderResult);
end

load(fullfile('kernels','Levin09.mat'));
kernelType = 1; % 1~8
if kernelType > 8
    k = fspecial('gaussian', 25, 1.6);
else
    k = kernels{kernelType};
end
sigmas      = [2, 2.55, 7.65]/255;
sigma       = sigmas(3);
totalIter   = 30; % default
lamda       = (sigma^2)/3; % default 3, ****** from {1 2 3 4} ******
modelSigma1 = 49; % default
modelSigma2 = 13; % ****** from {1 3 5 7 9 11 13 15} ******
modelSigmaS = logspace(log10(modelSigma1),log10(modelSigma2),totalIter);
rho         = sigma^2/((modelSigma1/255)^2);

ns          = min(25,max(ceil(modelSigmaS/2),1));
ns          = [ns(1)-1,ns];

for n_set = 1 : numel(setTest)
    %%% read images
    setTestCur = cell2mat(setTest(n_set));
    disp('--------------------------------------------');
    disp(['----',setTestCur,'-----Image Debluring-----']);
    disp('--------------------------------------------');
    folderTestCur = fullfile(folderTest,setTestCur);
    ext                 =  {'*.jpg','*.png','*.bmp'};
    filepaths           =  [];
    for i = 1 : length(ext)
        filepaths = cat(1,filepaths,dir(fullfile(folderTestCur, ext{i})));
    end
    eval(['PSNR_',taskTestCur,'_',setTestCur,' = zeros(length(filepaths),1);']);

    %%% folder to store results
    folderResultCur = fullfile(folderResult, ['Deblur_',setTestCur,'_kernel_',num2str(kernelType)]);
    if ~exist(folderResultCur,'file')
        mkdir(folderResultCur);
    end

    for i = 1 : length(filepaths)

        x  = imread(fullfile(folderTestCur,filepaths(i).name));
        [~,imageName,ext] = fileparts(filepaths(i).name);
        randn('seed',0);
        y = imfilter(im2double(x), k, 'circular', 'conv') + sigma*randn(size(x));
        [w,h,c]  = size(y);
        V = psf2otf(k,[w,h]);
        denominator = abs(V).^2;

        if c>1
            denominator = repmat(denominator,[1,1,c]);
            V = repmat(V,[1,1,c]);
        end
        upperleft   = conj(V).*fft2(y);

        if c==1
            load(fullfile(folderModel,'modelgray.mat'));
        elseif c==3
            load(fullfile(folderModel,'modelcolor.mat'));
        end
        z = single(y);
        if useGPU
            z           = gpuArray(z);
            upperleft   = gpuArray(upperleft);
            denominator = gpuArray(denominator);
        end
        tic;
        for itern = 1:totalIter
            %%% step 1
            rho = lamda*255^2/(modelSigmaS(itern)^2);
            z = real(ifft2((upperleft + rho*fft2(z))./(denominator + rho)));
            if ns(itern+1)~=ns(itern)
                [net] = loadmodel(modelSigmaS(itern),CNNdenoiser);
                net = vl_simplenn_tidy(net);
                if useGPU
                    net = vl_simplenn_move(net, 'gpu');
                end
            end
            %%% step 2
            res = vl_simplenn(net, z,[],[],'conserveMemory',true,'mode','test');
            residual = res(end).x;
            z = z - residual;
        end

        if useGPU
            output = im2uint8(gather(z));
        end
        toc;

        [PSNR_Cur,SSIM_Cur] = Cal_PSNRSSIM(x,output,0,0); %%% single
        disp(['Image Deblurring     ',num2str(PSNR_Cur,'%2.2f'),'dB','    ',filepaths(i).name]);
        eval(['PSNR_',taskTestCur,'_',setTestCur,'(',num2str(i),') = PSNR_Cur;']);

        if showResult
            imshow(cat(2,im2uint8(y),output,x));
            drawnow;
            title(['Image Deblurring     ',filepaths(i).name,'    ',num2str(PSNR_Cur,'%2.2f'),'dB'],'FontSize',12)
            pause(pauseTime)
            %pause()
            imwrite(output,fullfile(folderResultCur,[imageName,'.png']));
        end
    end
    disp(['Average PSNR is ',num2str(mean(eval(['PSNR_',taskTestCur,'_',setTestCur])),'%2.2f'),'dB']);

    %%% save PSNR
    save(fullfile(folderResultCur,['PSNR_',taskTestCur,'_',setTestCur,'.mat']),['PSNR_',taskTestCur,'_',setTestCur])

end

author: cszn

repo: IRCNN

save_path: github-repos/MATLAB/cszn-IRCNN


text:
function F = erfc(F, varargin)
%ERFC   Complementary error function of a CHEBFUN.
%   ERFC(X) is the complementary error function of the real-valued CHEBFUN X.
%   The complementary error function is defined as:
%       ERFC(X)(s) = 2/sqrt(pi) * integral from X(s) to inf of exp(-t^2) dt.
%                  = 1 - ERF(X)(s).
%
% See also ERF, ERFCX, ERFINV, ERFCINV.

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

% Input must be real:
if ( ~isreal(F) )
    error('CHEBFUN:CHEBFUN:erfc:notreal', 'Input must be real.');
end

% Call the compose method:
F = compose(F, @erfc, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = sampleinthist(varargin)
% VL_SAMPLINTHIST  Sample integral histogram
%   HISTS = VL_SAMPLINTHIST(INTHIST, BOXES) samples the integral
%   histogram INTHIST to obtain the histograms of the specified
%   BOXES.
%
%   INTHIST is a MxNxK array, where M x N are ``spatial'' dimensions,
%   and K is the number of histogram bins. INTHIST may be of class
%   UINT32 or DOUBLE.
%
%   Each box is a four dimensional vector [IMIN JMIN IMAX JMAX]' of
%   class UINT32 and correspond to the index set [IMIN, IMAX] x [JMIN,
%   JMAX]. To specify an empty box, let IMIN > IMAX.
%
%   HISTS stores one histogram per column (one for each box) and has K
%   rows, one for each histogram bin. HIST is of the same class of
%   INTHIST.
%
%   See also: VL_INTHIST(), VL_IMINTEGRAL(), VL_HELP().
[varargout{1:nargout}] = vl_sampleinthist(varargin{:});

author: yihui-he

repo: panorama

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


text:
% +TFUN
%
% Files

author: CelsoReyes


text:
function c = randchar(n)

u = rand(1,n);
u = floor(u*26)+65;
c=char(u);

author: bashtage

repo: mfe-toolbox


text:
% read_eeglabheader() - import EEGLAB dataset files
%
% Usage:
%   >> header = read_eeglabheader(filename);
%
% Inputs:
%   filename - [string] file name
%
% Outputs:
%   header   - FILEIO toolbox type structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2008-

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD, arno@salk.edu
%
% 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 header = read_eeglabheader(filename)

if nargin < 1
  help read_eeglabheader
  return
end

if ~isstruct(filename)
  s = whos('-file', filename);
  if any(strcmp({s.name}', 'EEG'))
    load('-mat', filename, 'EEG');
  else % EEGLAB > 2021.0 saves content of EEG as default
    EEG = load('-mat', filename);
  end
else
  EEG = filename;
end

header.Fs          =  EEG.srate;
header.nChans      =  EEG.nbchan;
header.nSamples    =  EEG.pnts;
header.nSamplesPre = -EEG.xmin*EEG.srate;
header.nTrials     =  EEG.trials;
try
  header.label     = { EEG.chanlocs.labels }';
catch
  ft_warning('creating default channel names');
  for i=1:header.nChans
    header.label{i} = sprintf('chan%03d', i);
  end
end
ind = 1;
for i = 1:length( EEG.chanlocs )
  if isfield(EEG.chanlocs(i), 'X') && ~isempty(EEG.chanlocs(i).X)
    header.elec.label{ind, 1} = EEG.chanlocs(i).labels;
    % this channel has a position
    header.elec.elecpos(ind,1) = EEG.chanlocs(i).X;
    header.elec.elecpos(ind,2) = EEG.chanlocs(i).Y;
    header.elec.elecpos(ind,3) = EEG.chanlocs(i).Z;
    ind = ind+1;
  end
end

if ind>1 && isfield(EEG, 'chaninfo') && isfield(EEG.chaninfo, 'nosedir') && ~isempty(EEG.chaninfo.nosedir)
  switch EEG.chaninfo.nosedir
    case '+X'
      header.elec.coordsys = 'ctf';
    case '-X'
      header.elec.coordsys = 'pls';
    case '+Y'
      header.elec.coordsys = 'ras';
    case '-Y'
      header.elec.coordsys = 'lpi';
    otherwise
  end
end

% remove data
% -----------
%if isfield(EEG, 'datfile')
%    if ~isempty(EEG.datfile)
%        EEG.data = EEG.datfile;
%    end
%else
%    EEG.data = 'in set file';
%end;
EEG.icaact = [];

header.orig = EEG;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [signalPeaks, signalPeaksArray, signalSigmas, signalStruct] = computeSignalPeaks(signalMatrix, varargin)
	% [signalPeaks, signalPeaksArray, signalSigmas] = computeSignalPeaks(signalMatrix, varargin)
	%
	% Binarize [0,1] input analog signals based on peaks in the signal.
	%
	% Biafra Ahanonu
	% started: 2013.10.28
	%
	% inputs
	% 	signalMatrix: [nSignals frame] matrix containing analog input signals.
	% outputs
	% 	signalPeaks: [nSignals frame] matrix. Binary matrix with 1 = peaks, 0 = non-peaks.
	% 	signalPeaksArray: {1 nSignals} cell array. Each cell contains [1 nPeaks] vector that stores the frame locations of each peak.
	% 	signalSigmas: [nSignals 1] - std of each signal.
	% 	signalStruct: structure containing signalPeaks and signalPeaksArray if multiple thresholds requested.
	% options
	% 	See below.
	% 	% make a plot?
	% 	options.makePlots = 0;
	% 	% show waitbar?
	% 	options.waitbarOn = 1;
	% 	% make summary plots of spike information
	% 	options.makeSummaryPlots = 0;
	% 	% number of standard deviations above the threshold to count as spike
	% 	options.numStdsForThresh = 3;
	% 	% minimum number of time units between events
	% 	options.minTimeBtEvents = 8;
	% 	% shift peak detection
	% 	options.nFramesShift = 0;
	% 	% should diff and fast oopsi be done?
	% 	options.addedAnalysis = 0;
	% 	% use simulated oopsi data
	% 	options.oopsiSimulated = 0;
	% changelog
		% 2015.10.06 [00:14:09] Changed computePeakForSignal to shift the signal to the actual nearby peak since findpeak is sometimes off by a frame or two, should also improve the S-ratio.
		% 2016.07.05 [14:52:43] Made changes to computePeakForSignal to improve diff based peak detection.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.04.20 [21:29:09] - Better comments for options.
	% TODO:
		% Add option to obtain multiple signal peak outputs from different thresholds in the same run, would save time.
		% allow input of options file (e.g. for different GCaMP variants, brain regions, etc.)
		% integrate nearest neighbor into analysis if there is a lot of cross-talk
		% possibly integrate into identifySpikes?
		% TODO: convert main loop to parfor, convert signalMatrix to cell array to allow this. Normally fast enough that this won't provide a speed-up.

	% add controller directory and subdirectories to path
	% addpath(genpath(pwd));
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% Binary: 1 = show plots with found events and other information for each signal. 0 = do not show signal plot GUI.
	options.makePlots = 0;
	% Binary: 1 = show wait bar, 0 = do not show wait bar.
	options.waitbarOn = 1;
	% Binary: 1 = make summary plots of spike information, 0 = no summary plots.
	options.makeSummaryPlots = 0;
	% ===
	% Int: number of standard deviations above the threshold to count as spike
	options.numStdsForThresh = 3; % 0.5
	% DEPRECIATED - Int: alternative to options.numStdsForThresh for display purposes
	options.numStdsForThreshTwo = 2;
	% Int: minimum number of time units between events.
	options.minTimeBtEvents = 8;
	% Str: detect on differential ('diff') or raw ('raw') trace
	options.detectMethod = 'diff'; % 'raw'
	% Int: the size of the window to use to ignore smaller peaks near larger peaks.
	options.movAvgReqSize = 2;
	% Int: the size of the moving average to use on the input signals, to smooth out noise.
	options.movAvgFiltSize = 3;
	% Binary: 1 = use moving average as specified in options.movAvgFiltSize.
	options.doMovAvg = 1;
	% Binary: 1 = subtract median calculated over a filter of some range.
	options.doMedianFilter = 1;
	% Int: number of frames to calculate rolling median filter.
	options.medianFilterLength = 201;
	% Binary: 1 = leave report peaks as determined by findpeaks (e.g. normally midpoint of the peak rise if using 'diff'). 0 = adjust peak location to the maximum value found within a options.peakMaxLook window.
	options.reportMidpoint = 0;
	% Int vector: frames before and after region around each peak to look for a maximum to adjust the test peak by.
	options.peakMaxLook = -6:6;
	% Int: number of frames to shift detected peaks.
	options.nFramesShift = 0;
	% ===
	% Binary: 1 = perform diff and fast oopsi.
	options.addedAnalysis = 0;
	% Binary: 1 = use simulated oopsi data.
	options.oopsiSimulated = 0;
	% Binary: 1 = open workers, 0 = do not open workers
	options.parallel = 1;
	% Binary: 1 = display output information. 0 = suppress most output (e.g. when running in batch for certain applications).
	options.outputInfo = 1;
	% Binary: 1 = convert input inputSignals matrix to cell array
	options.convertSignalsToCell = 1;
	% get user inputs
	options = getOptions(options,varargin);
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	%     eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================
	try
		% make sure matrix is double
		% if ~strcmp(class(signalMatrix),'double')
		if ~isa(signalMatrix,'double')
			if options.outputInfo==1
				disp('converting signalMatrix to double');
			end
			signalMatrix = double(signalMatrix);
		end
		if options.outputInfo==1
			disp('calculating signal peaks...')
		end
		nSignals = size(signalMatrix,1);
		nFrames = size(signalMatrix,2);

		% Remove any NaN signals, set to the mean
		if sum(isnan(signalMatrix(:)))>0
			if options.outputInfo==1
				disp('Removing NaNs from matrix');
			end
			for signalNo = 1:nSignals
				repIdx = isnan(signalMatrix(signalNo,:));
				signalMatrix(signalNo,repIdx) = nanmean(signalMatrix(signalNo,:));
			end
		end

		signalStruct = struct;

		% this matrix will contain binarized version of signalMatrix
		signalPeaks = zeros(size(signalMatrix));
		% contains a list for each signal of locations of peaks
		signalPeaksArray = cell([1 nSignals]);
		% open waitbar
		% if options.waitbarOn==1 waitbarHandle = waitbar(0, 'detecting traces...'); end
		% loop over all signals.
		% reverseStr = '';
		manageParallelWorkers('parallel',options.parallel);
		signalPeaksArrayTmp = cell([1 nSignals]);
		signalSigmas = repmat([],[1 nSignals]);
		% try;[percent progress] = parfor_progress(nSignals);catch;end; dispStepSize = round(nSignals/20); dispstat('','init');

		% Convert to cell array to reduce memory transfer during parallelization
		if options.convertSignalsToCell==1
			% signalMatrix = squeeze(mat2cell(signalMatrix,ones(1,size(signalMatrix,1)),size(signalMatrix,2)));
		end

		[optionsOut] = computePeakForSignalOptions('options', options);

		% Only implement in Matlab 2017a and above
		if ~verLessThan('matlab', '9.2')
			D = parallel.pool.DataQueue;
			afterEach(D, @nUpdateParforProgress);
			p = 1;
			N = nSignals;
			nInterval = round(nSignals/30);%25
			options_waitbarOn = options.waitbarOn;
		end

		% signalIdxAll = {};
		% if options.convertSignalsToCell==1
		%     parfor signalNum3 = 1:nSignals
		%         signalIdxAll{signalNum3} = signalNum3;
		%     end
		% else
		%     parfor signalNum3 = 1:nSignals
		%         signalIdxAll{signalNum3} = sub2ind(size(signalMatrix), repmat(signalNum3,[nFrames 1]), [1:nFrames]');
		%     end
		% end

		% startState = ticBytes(gcp);
		if isempty(gcp)
			optsConstant.Value = optionsOut;
		else
			optsConstant = parallel.pool.Constant(optionsOut);
		end
		optionsCopy_addedAnalysis = options.addedAnalysis;
		optionsCopy_convertSignalsToCell = options.convertSignalsToCell;
		parfor signalNum = 1:nSignals
			optionsOutCopy = optsConstant.Value;
			% [percent progress] = parfor_progress;if mod(progress,dispStepSize) == 0;dispstat(sprintf('progress %0.1f %',percent));else;end
			% get the current signal and find its peaks
			% if options.convertSignalsToCell==1
			%     thisSignal = signalMatrix{signalNum2};
			%     signalNum = signalNum2;
			%     signalNum3 = signalNum2;
			%     signalSigmas(signalNum2) = std(thisSignal);
			%     signalPeaksArray{signalNum2} = computePeakForSignal(thisSignal,optionsOutCopy);
			% else
			%     signalNum = find(cellfun(@(x) sum(signalNum2==x),signalIdxAll));
			%     % signalIdx = sub2ind(size(signalMatrix), repmat(signalNum,[nFrames 1]), [1:nFrames]');
			%     % signalIdx = signalIdxAll{signalNum};
			%     thisSignal = signalMatrix(signalNum2);
			%     signalSigmas(signalNum) = std(thisSignal);
			%     signalPeaksArray{signalNum} = computePeakForSignal(thisSignal,optionsOutCopy);
			% end

			thisSignal = signalMatrix(signalNum,:);

			% if iscell(signalMatrix)==1
			%     thisSignal = signalMatrix{signalNum};
			% else
			%     thisSignal = signalMatrix(signalNum,:);
			%     % signalSigmas(signalNum) = std(thisSignal);
			%     % signalPeaksArray{signalNum} = computePeakForSignal(thisSignal,optionsOutCopy);
			% end

			%
			signalSigmas(signalNum) = std(thisSignal);
			signalPeaksArray{signalNum} = computePeakForSignal(thisSignal,optionsOutCopy);

			% [~] = viewComputePeaksPlot(thisSignal,signalPeaksArray{signalNum},[0 0 0],options.makePlots,50,2,'on')
			% ===
			if optionsCopy_addedAnalysis==1
				% using diff
				detectOld = optionsOutCopy.detectMethod;
				optionsOutCopy.detectMethod = 'diff';
				optionsOutCopy.detectMethod = detectOld;

				signalPeaksArrayTmp{signalNum} = computePeakForSignal(thisSignal,optionsOutCopy);
				% plot the resulting peaks overlayed
				[~] = viewComputePeaksPlot(thisSignal,signalPeaksArrayTmp{signalNum},[0 0 1],options.makePlots,20,0.1,'off',options)
				legend('signal','raw','signal','diff')
				title(['raw: ' num2str(length(signalPeaksArray{signalNum})) ' | diff: ' num2str(length(signalPeaksArrayTmp{signalNum}))]);
				[~,~,~]=ginput(1);
				% fast oopsi
				[Nhat] = computePeakForSignalOopsi(thisSignal,signalPeaksArray{signalNum},options);
				[~,~,~]=ginput(1);
				% if options.convertSignalsToCell==1
				% else
				% end
				% ===
			end
			% create matrix of peaks
			% signalPeaks(signalNum,signalPeaksArray{signalNum})=1;

			% waitbar access
			% reverseStr = cmdWaitbar(signalNum,nSignals,reverseStr,'inputStr','detecting peaks','waitbarOn',options.waitbarOn,'displayEvery',50);

			if ~verLessThan('matlab', '9.2')
				% Update
				send(D, signalNum);
			end
		end
		% tocBytes(gcp,startState)
		if options.makePlots==1
			for signalNum=1:nSignals
				if options.convertSignalsToCell==1
					% thisSignal = signalMatrix{signalNum};
					thisSignal = signalMatrix(signalNum,:);
				else
					thisSignal = signalMatrix(signalNum,:);
				end
					[~] = viewComputePeaksPlot(thisSignal,signalPeaksArray{signalNum},[0 0 0],options.makePlots,50,2,'on',options,signalNum,nSignals,options.numStdsForThresh);
					pause
					keyIn = get(gcf,'CurrentCharacter');
					if strcmp(keyIn,'e')
						return
					end
				% pause
				% clf
			end
		end
		for signalNum=1:nSignals
			% create matrix of peaks
			signalPeaks(signalNum,signalPeaksArray{signalNum})=1;
		end
	catch err
		signalPeaks = [];
		signalPeaksArray = {};
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
	% summary of general statistics for this set of IC data
	try
		if options.makeSummaryPlots==1
			viewSpikeSummary(signalMatrix,signalPeaks);
		end
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
	function nUpdateParforProgress(~)
		if ~verLessThan('matlab', '9.2')
			p = p + 1;
			if (mod(p,nInterval)==0||p==2||p==nSignals)&&options_waitbarOn==1
				if p==nSignals
					fprintf('%d\n',round(p/nSignals*100))
				else
					fprintf('%d|',round(p/nSignals*100))
				end
				% cmdWaitbar(p,nSignals,'','inputStr','','waitbarOn',1);
			end
			% [p mod(p,nInterval)==0 (mod(p,nInterval)==0||p==nSignals)&&options_waitbarOn==1]
		end
	end
end
function [inputSignal] = viewComputePeaksPlot(inputSignal,testpeaks,dotColor,makePlots,markersize,linewidth,holdVal,options,signalNum,nSignals,numStdsForThresh)
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% decide whether to plot the peaks with points indicating location of
	% chosen peaks
	if makePlots==1
		% setFigureDefaults()
		fig1 = figure(422); clf;
		subplot(4,3,1)
			% hist(inputSignal(testpeaks),20);
			hist((inputSignal(testpeaks) - nanmean(inputSignal(:)))/nanstd(inputSignal(:)),20);
			hold(holdVal);
			title(num2str(nanstd(inputSignal(:))))
			xlabel('Peak amplitude (Z-score)');ylabel('Peak count')
			set(gca,'XMinorTick','on','TickDir','out');box off;
			% plot(histBins,histCounts);
			% set(gca,'yscale','log');

		subplot(4,3,2)
			peakROI = [-20:20];
			extractMatrix = bsxfun(@plus,testpeaks',peakROI);
			extractMatrix(extractMatrix<=0)=1;
			extractMatrix(extractMatrix>=size(inputSignal,2))=size(inputSignal,2);
			spikeCenterTrace = reshape(inputSignal(extractMatrix),size(extractMatrix));
			plot(repmat(peakROI, [size(spikeCenterTrace,1) 1])', spikeCenterTrace','Color',[4 4 4]/8)
			set(gca,'TickDir','out');box off;
			xlabel('Time (frames)');ylabel('Signal amplitude')
			title('All peaks')

		subplot(4,3,3)
			peakSignalAmplitude = inputSignal(testpeaks(:));
			[peakSignalAmplitude, peakIdx] = sort(spikeCenterTrace(:,round(end/2)+1),'descend');
			spikeCenterTrace = spikeCenterTrace(peakIdx,:);
			if size(spikeCenterTrace,1)>20
				spikeCenterTrace = spikeCenterTrace(1:20,:);
			end
			% xlabel('Time (frames)')
			plot(repmat(peakROI, [size(spikeCenterTrace,1) 1])', spikeCenterTrace','Color',[4 4 4]/8)
			xlabel('Time (frames)');ylabel('Signal amplitude')
			title('Top 20 peaks')
			set(gca,'TickDir','out');box off;

		subplot(4,3,[4:12])
			set(gcf,'color','w');
			% scnsize = get(0,'ScreenSize');
			% position = get(fig1,'Position');
			% outerpos = get(fig1,'OuterPosition');
			% borders = outerpos - position;
			% edge = -borders(1)/2;
			% %pos1 = [scnsize(3)/2 + edge, 0, scnsize(3)/2 - edge, scnsize(4)];
			% pos1 = [0, 0, scnsize(3), scnsize(4)];
			% set(fig1,'OuterPosition',pos1);

			plot(inputSignal, 'r');

			inputSignalMedian=medfilt1(inputSignal,options.medianFilterLength,'omitnan','truncate');
			inputSignal = inputSignal - inputSignalMedian;
			inputSignal = filtfilt(ones(1,options.movAvgFiltSize)/options.movAvgFiltSize,1,inputSignal);
			hold on;
			plot(inputSignal, 'b');

			set(gca,'Color','none'); box off;
			set(gca,'TickDir','out');box off;
			hold on;
			% axis([0 length(inputSignal) -0.1 0.5]);
			scatter(testpeaks, inputSignal(testpeaks),markersize, 'LineWidth',linewidth,'MarkerFaceColor',dotColor, 'MarkerEdgeColor',dotColor)
			legend({'Raw signal','Raw signal minus rolling median filter'},'Location','northoutside')
			axis tight;
			xlabel('Time (frames)');ylabel('Signal amplitude')
			% options.numStdsForThresh = options.numStdsForThreshTwo;
			% signalPeaksArray{signalNum} = computePeakForSignal(thisSignal, 'options', options);

			% [x,y,reply]=ginput(1);
			% close(fig1);
			hold(holdVal);
			title(sprintf('%d peaks | Zoom is enabled for closer look at peaks',length(testpeaks)))
		zoom on
		ciapkg.overloaded.suptitle(sprintf('Right arrow key to move to next signal | press ''e'' to exit | threshold = %0.2f | signal %d/%d',numStdsForThresh,signalNum,nSignals))
	end
end
function [Nhat] = computePeakForSignalOopsi(inputSignal,testpeaks, options, varargin)
	% clear, clc,
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	switch options.oopsiSimulated
		case 1
			% set simulation metadata
			T       = 1000; % # of time steps
			V.dt    = 1/8;  % time step size

			% initialize params
			P.a     = 1;    % observation scale
			P.b     = 0;    % observation bias
			tau     = 1.5;    % decay time constant
			P.gam   = 1-V.dt/tau; % C(t) = gam*C(t-1)
			P.lam   = 0.1;  % firing rate = lam*dt
			P.sig   = 0.1;  % standard deviation of observation noise

			% simulate data
			N = poissrnd(P.lam*V.dt*ones(T,1)); % simulate spike train
			N(100) = 10;
			P.sig   = std(N);  % standard deviation of observation noise
		case 0
			% our data
			% inputSignal = inputSignal(1:1000);
			% testpeaks = testpeaks(testpeaks<1000);
			T = length(inputSignal);
			V.dt    = 1/5;  % time step size
			P.k = 30;
			%
			P.a     = 1;    % observation scale
			P.b     = 0;    % observation bias
			tau     = 1.5;    % decay time constant
			P.gam   = 1-V.dt/tau; % C(t) = gam*C(t-1)
			P.lam   = 0.01;  % firing rate = lam*dt

			inputSignalMetric = nanmean(inputSignal(testpeaks));
			% inputSignalMetric = nanmax(inputSignal(testpeaks));
			N = inputSignal(:)/inputSignalMetric;
			% inputSignal = normalizeVector(inputSignal,'normRange','zeroToOne');
			% N = inputSignal(:);
			P.sig = std(inputSignal(:)); % standard deviation of observation noise
			P.sig
			F = N;
			% N = zeros([length(inputSignal) 1]);
			% N(testpeaks) = 1;
			% normalizeVector(inputSignal(:),'normRange','zero')

		otherwise
			% body
	end
	P.sig = 0.1; % standard deviation of observation noise
	C = filter(1,[1 -P.gam],N);         % calcium concentration
	F = P.a*C+P.b + P.sig*randn(T,1);   % observations
	% fast oopsi
	[Nhat Phat] = fast_oopsi(F,V,P);

	% smc-oopsi
	V.smc_iter_max = 1;
	% [M P V] = smc_oopsi(F,V,P);

	%% plot results
	figure(1), clf
	tvec=0:V.dt:(T-1)*V.dt;
	tvec = 1:T;
	h(1)=subplot(411); plot(tvec,F); axis('tight'), ylabel('F (au)')
	h(2)=subplot(412); plot(tvec,C); axis('tight'), ylabel('C (au)')
	switch options.oopsiSimulated
		case 1
			h(3)=subplot(4,1,[3 4]); stem(tvec,N,'.'); hold on, plot(tvec,Nhat,'r','linewidth',1), axis('tight'), ylabel('fast')
		case 0
			h(3)=subplot(4,1,3); stem(tvec,N,'.'); hold on, plot(tvec,Nhat,'r','linewidth',1), axis('tight'), ylabel('fast')
			h(4)=subplot(414);plot(tvec,inputSignal, 'r'); box off; hold on;
			scatter(testpeaks, inputSignal(testpeaks),50, 'LineWidth',2,'MarkerFaceColor',[0 0 1], 'MarkerEdgeColor',[0 0 1]); axis('tight')
		otherwise
	end
	% Nsmc = M.nbar/max(M.nbar);
	% Nsmc(Nsmc<0.1)=0;
	% h(4)=subplot(414); stem(tvec,N); hold on, plot(tvec,Nsmc,'k','linewidth',2); axis('tight'), ylabel('smc')
	% xlabel('time (sec)')
	linkaxes(h,'x')
end
function [optionsOut] = computePeakForSignalOptions(varargin)
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% number of standard deviations above the threshold to count as spike
	options.numStdsForThresh = 3;
	% minimum number of time units between events
	options.minTimeBtEvents = 10;
	% make a plot?
	options.makePlots = 0;
	% the size of the moving average
	options.movAvgReqSize = 2;
	options.movAvgFiltSize = 3;
	% subtract median calculated over a filter of some range.
	options.doMedianFilter = 1;
	% number of frames to calculate median filter
	options.medianFilterLength = 201;
	% decide whether to have a moving average
	options.doMovAvg = 1;
	% report the midpoint of the rise
	options.reportMidpoint=0;
	% shift peak detection
	options.nFramesShift = 0;
	% detect on differential ('diff') or raw ('raw') trace
	options.detectMethod = 'raw';
	% region around each peak to look for a maximum to adjust the test peak by
	options.peakMaxLook = [-6:6];
	% get options
	optionsOut = getOptions(options,varargin,'showWarnings',0);
end
function [testpeaks] = computePeakForSignal(inputSignal, options)
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% identifies peaks in an input signal given a particular threshold and other parameters.
	% biafra ahanonu
	% started: 2013.10.28
	% adapted from Lacey Kitch's and Laurie Burns' code
	% inputs
		%
	% outputs
		%
	% changelog
		% 2013.11.18 [20:06:00]
	% TODO
		%

	%========================
	% unpack options into current workspace
	fn=fieldnames(options);
	for i=1:length(fn)
		eval([fn{i} '=options.' fn{i} ';']);
	end;
	%========================

	% histBins = logspace(1,max(inputSignal));
	% [histCounts histBins] = hist(inputSignal(:),100);
	% plot(histBins,histCounts);
	% set(gca,'yscale','log');

	% median filter subtract
	if doMedianFilter==1
		inputSignalMedian=medfilt1(inputSignal,options.medianFilterLength,'omitnan','truncate');
		inputSignal = inputSignal - inputSignalMedian;
	end
	% moving average of input signal
	if doMovAvg==1
		% class(movAvgFiltSize)
		% movAvgFiltSize
		% class(inputSignal)
		% inputSignal(1:10)
		inputSignal = filtfilt(ones(1,movAvgFiltSize)/movAvgFiltSize,1,inputSignal);
	end

	switch options.detectMethod
		case 'diff'
			rawInputSignalStd = std(inputSignal(:));
			rawInputSignal = inputSignal;
			% get the differential
			inputSignal = [0 diff(inputSignal)];
			inputSignal(inputSignal<0) = 0;
			% options.nFramesShift = 0;
		case 'raw'
			%
		otherwise
			% body
	end

	% get standard deviation of current signal
	% inputSignalStd = std(inputSignal(:));
	inputSignalStd = nanstd(inputSignal(:));
	thisStdThreshold = inputSignalStd*numStdsForThresh;

	% =======
	% peakIdx = bsxfun(@plus,options.timeSeq',testpeaks);
	% tmpTestPeak = testpeaks;
	% % remove peaks outside range of signal
	% peakIdx(peakIdx>length(loopSignal))=[];
	% peakIdx(peakIdx<=0)=[];

	% % remove signal then add back in noise based on signal statistics
	% noiseSignal = loopSignal;
	% noiseSignal(peakIdx) = NaN;
	% =======

	% run findpeaks (part of signal), returns maxima above thisStdThreshold
	% and ignores smaller peaks around larger maxima within minTimeBtEvents
	warning off
	[~,testpeaks] = findpeaks(inputSignal,'minpeakheight',thisStdThreshold,'minpeakdistance',minTimeBtEvents);
	warning on

	% ignores smaller peaks around larger maxima within minTimeBtEvents
	%[~,testpeaks2] = findpeaks(inputSignal,'minpeakdistance',minTimeBtEvents);
	%testpeaks = intersect(testpeaks,testpeaks2);
	% extra check
	testpeaks = intersect(testpeaks,...
		find(filtfilt(ones(1,movAvgReqSize)/movAvgReqSize,1,inputSignal)>thisStdThreshold)...
		);

	% remove non-peaks within some set criteria
	switch options.detectMethod
		case 'diff'
			% get the differential
			testpeaks = testpeaks(rawInputSignal(testpeaks)>(numStdsForThresh*rawInputSignalStd));
			% options.nFramesShift = 0;

			% switch back for later analysis
			inputSignal = rawInputSignal;
		case 'raw'
			%
		otherwise
			% body
	end

	% check that maximum is at peak, else shift it
	if ~isempty(testpeaks)&options.reportMidpoint==0
		peakMaxLook = options.peakMaxLook;
		for peakNo = 1:length(testpeaks)
				testNewPeakIdx = testpeaks(peakNo)+peakMaxLook;
				if min(testNewPeakIdx)<=0|max(testNewPeakIdx)>length(inputSignal)
					continue;
				end
				testSignal = inputSignal(testNewPeakIdx);
				[~,maxSignal] = max(testSignal);
				% [maxSignal testpeaks(peakNo) testpeaks(peakNo)+(maxSignal-round(length(testNewPeakIdx)/2))]
				testpeaks(peakNo) = testpeaks(peakNo)+(maxSignal(1)-round(length(testNewPeakIdx)/2));
		end
	elseif options.reportMidpoint==1
		% peakMaxLook = options.peakMaxLook;
		% for peakNo = 1:length(testpeaks)
		%         testNewPeakIdx = testpeaks(peakNo)+peakMaxLook;
		%         if min(testNewPeakIdx)<=0|max(testNewPeakIdx)>length(inputSignal)
		%             continue;
		%         end
		%         testSignal = inputSignal(testNewPeakIdx);
		%         [~,maxSignal] = max(testSignal);
		%         % [maxSignal testpeaks(peakNo) testpeaks(peakNo)+(maxSignal-round(length(testNewPeakIdx)/2))]
		%         testpeaks(peakNo) = testpeaks(peakNo)+(maxSignal(1)-round(length(testNewPeakIdx)/2));
		% end
	end

	% shift peaks
	testpeaks = testpeaks + options.nFramesShift;

end

function [peakIdx] = subfxnCalcSignalNew(noiseSigmaThreshold,noiseSignal,loopSignal,tmpTestPeak)
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	noiseStd = nanstd(noiseSignal(:));
	% noiseStd
	nPeaks = length(tmpTestPeak);
	peakIdxTmp1 = cell([1 nPeaks]);
	peakIdxTmp2 = cell([1 nPeaks]);
	% noiseSigmaThreshold = options.noiseSigmaThreshold;
	loopSignalThresholded = loopSignal>noiseSigmaThreshold*noiseStd;
	loopSignalThresholdedDiff = [0 diff(loopSignalThresholded)];
	for peakNo = 1:nPeaks
		try
			peakFrame = tmpTestPeak(peakNo);

			% currFrame = peakFrame;
			signalCriteria = loopSignalThresholdedDiff(1:peakFrame);
			currFrame = find(signalCriteria,1,'last');
			if isempty(currFrame)
				if sum(loopSignalThresholded(1:peakFrame))==length(loopSignalThresholded(1:peakFrame))
					currFrame = 1;
				end
			end
			% currFrame
			% aboveNoise = 1;
			% while aboveNoise==1
			%   aboveNoise = loopSignal(currFrame)>options.noiseSigmaThreshold*noiseStd;
			%   currFrame = currFrame - 1;
			%   if currFrame<1
			%       break
			%   end
			% end
			peakIdxTmp1{peakNo} = currFrame:peakFrame;

			% currFrame = peakFrame;
			% loopSignalThresholded = loopSignal>noiseSigmaThreshold*noiseStd;
			signalCriteria = abs(loopSignalThresholdedDiff(peakFrame:end));
			currFrame = find(signalCriteria,1,'first')+peakFrame;
			if isempty(currFrame)
				if sum(loopSignalThresholded(peakFrame:end))==length(loopSignalThresholded(peakFrame:end))
					currFrame = length(loopSignal);
				end
			end
			% currFrame
			% aboveNoise = 1;
			% currFrame = peakFrame;
			% while aboveNoise==1
			%   aboveNoise = loopSignal(currFrame)>options.noiseSigmaThreshold*noiseStd;
			%   currFrame = currFrame + 1;
			%   if currFrame>length(loopSignal)
			%       break
			%   end
			% end
			peakIdxTmp2{peakNo} = peakFrame:currFrame;
		catch err
			fprintf('peakFrame = %d\n',peakFrame);
			display(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			display(repmat('@',1,7))
		end
	end
	peakIdx = [peakIdxTmp1{:} peakIdxTmp2{:}];
	peakIdx = unique(peakIdx(:));
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function S = AccSign(p)
%ACCSIGN      Computes sign(sum(p))
%
%   S = AccSign(p)
%
%On return, S is the sign of sum(p), also in the presence
%  of underflow. Input vector p may be single or double precision.
%
%Implements Algorithm 8.2 from
%  S.M. Rump, T. Ogita, S. Oishi: Accurate Floating-point Summation II:
%    Sign, K-fold Faithful and Rounding to Nearest, Siam J. Sci. Comput.,
%    31(2):1269-1302, 2008.
%Requires (4m+2)n flops for m executions of repeat-until loop in Transform.
%
%Reference implementation! Slow due to interpretation!
%

% written  03/03/07     S.M. Rump
% modified 05/09/09     S.M. Rump  rounding to nearest, complex input
%

  if ~isreal(p)
    error('AccSign for real input only')
  end

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

  if isa(p,'double')
    nmax = 2^52;            % nmax = 450,359,962,737,0496
  else
    nmax = 2^23;            % nmax = 8,388,608
  end
  if length(p)>nmax
    error(['maximum length of input vector for AccSign ' int2str(nmax) '.'])
  end

  kPhi = 1;
  [tau1,tau2,p] = Transform(p,0,kPhi);
  S = sign(tau1);

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function SimpleWHPlotSpectrograms(W,H,Data, plotAll);
% plots seqNMF factors W and H
% Also plots Data if provided, and reconstruction if data is not provided
% plotAll=1 means plot all data
% plotAll=0 means crop data so it's on the same timescale as W's
% Emily Mackevicius 2/3/2018

clf

% set(gcf, 'color', 'w');
if nargin < 4
    plotAll = 1;
end
if nargin < 3 || size(Data,2)==0
    plotData = 0;
else
    plotData = 1;
end

[N,K,L] = size(W);
[~,T] = size(H);
color_palet = [[0 .6 .3]; [.7 0 .7]; [1 .6 0];  [.1 .3 .9];  [1 .1 .1];  [0 .9 .3]; [.4 .2 .7]; [.7 .2 .1]; [.1 .8 1 ]; [1 .3 .7]; [.2 .8 .2]; [.7 .4 1]; [.9 .6 .4]; [0 .6 1]; [1 .1 .3]];
% if K<7
%     color_palet = jet(K);
% else
%     color_palet = [jet(7); (jet(7)+1)/2];
%     color_palet = repmat(color_palet, ceil(K/size(color_palet,1)),1);
% end
kColors = color_palet(1:K,:);
%% set widths of subplots
m = .1; % margin
ww = .35; % width of W plot
wwflat = .05; % width of Wflat plot
hh = .45; % height of H plot
hdata = 1-hh-2*m;
wdata = 1-ww-wwflat-2*m;
sep = ceil(L*.1);

%% crop data, unless plotAll
if plotAll
    indplot = 1:T;
else
    indplot =(1:ceil((K*(L+sep))/ww*wdata)); % so that data and W's are on same scale
    indplot(indplot>T) = [];
end

%% plot W's
axW = subplot('Position', [m m ww hdata]);
hold on
set(gca, 'ColorOrder', kColors);

WsToPlot = zeros(N,K*(L+sep));
XsToPlot = zeros(3,K);
YsToPlot = [N*ones(1,K); zeros(2,K)]+.5;
for ki = 1:K
    WsToPlot(:,((L+sep)*(ki-1)+1):((L+sep)*(ki-1)+L)) = squeeze(W(:,ki,:));
    XsToPlot(:,ki) = [(L+sep)*(ki-1)+1 (L+sep)*(ki-1)+1 (L+sep)*(ki-1)+L];
end
clims = [0 prctile(WsToPlot(WsToPlot>0),99)]; % if all W's are empty this line will bug
imagesc(WsToPlot, clims);
cmap = 1/256*flipud([158,1,66;213,62,79;244,109,67;253,174,97;254,224,139;255,255,191;230,245,152;171,221,164;102,194,165;50,136,189;94,79,162;1 1 1]);
% cmap = flipud(gray);
colormap(cmap)
plot(XsToPlot,YsToPlot, 'linewidth', 2);
xlim([1 K*(L+sep)]);ylim([0 N+.1]+.5)
set(gca, 'ydir', 'reverse')
axis off

%% plot data
axIm = subplot('Position', [m+ww m wdata hdata]);
if plotData
    clims = [0 prctile(Data(Data>0),99)];
    Im = imagesc(Data(:,indplot), clims);
else
    toplot = helper.reconstruct(W,H(:,indplot));
    clims = [0 prctile(toplot(:),99.9)];
    Im = imagesc(toplot,clims);
end
set(gca,'ydir','reverse')
hold on; plot([0 0 length(indplot)+1], [N 0 0]+.5, 'k')
xlim([0 length(indplot)+1]);ylim([0 N+.1]+.5)
axis off
%% plot Wflat (collapse out L dimension of W)
axWflat = subplot('Position', [m+ww+wdata m wwflat hdata]);
hold on
set(gca, 'ColorOrder', kColors);
plot(squeeze(sum(W,3)), 1:N,'>', 'markersize', 2.5);
ylim([0 N+.1]+.5)
axis tight
xlims = xlim;
xlim([xlims(2)*.1 xlims(2)])
set(gca, 'ydir', 'reverse')
axis off
%% plot H's
axH = subplot('Position', [m+ww m+hdata wdata hh]);
Hrescaled = squeeze(max(sum(W,1),[],3))'.*H; % rescale by approximate loading
dn = prctile(Hrescaled(:),100)/2;
for ki = K:-1:1
    Xs = [1 1:length(indplot) length(indplot)];
    Ys = [dn*ki (dn*ki + Hrescaled(K-ki+1,indplot)) dn*ki]-dn/2;
    patch(Xs,Ys, kColors(K-ki+1,:), 'edgecolor', kColors(K-ki+1,:))
    hold on
end
ylim([0 dn*K+dn*3]);xlim([0 length(indplot)+1])
axis off
%%
if plotAll
      linkaxes([axIm axW axWflat], 'y'); linkaxes([axIm axH], 'x');
end

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
function map = tab20c(N)
% Qualitative colormap from MatPlotLib, for plots using the line ColorOrder.
% In MatPlotLib 2 it is named VEGA20C, for MatPlotLib 3 was renamed TAB20C.
%
% Copyright (c) 2017-2020 Stephen Cobeldick
%
%%% Syntax:
%  map = tab20c
%  map = tab20c(N)
%
% For MatPlotLib 2.0 improved colormaps were created for plot lines of
% categorical data. The new colormaps are introduced here:
% <http://matplotlib.org/2.0.0rc2/users/dflt_style_changes.html>
% VEGA10/TAB10 is the default Line Color Order for MatPlotLib 2 and 3.
%
% MATLAB axes ColorOrder (note that this is NOT the axes COLORMAP):
% <https://www.mathworks.com/help/matlab/creating_plots/defining-the-color-of-lines-for-plotting.html>
% <https://www.mathworks.com/help/matlab/graphics_transition/why-are-plot-lines-different-colors.html>
%
%% Examples %%
%
%%% PLOT using matrices:
% N = 20;
% axes('ColorOrder',tab20c(N),'NextPlot','replacechildren')
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*sin(x+2*n*pi/N), X(:), 1:N);
% plot(X,Y, 'linewidth',4)
%
%%% PLOT in a loop:
% N = 20;
% set(0,'DefaultAxesColorOrder',tab20c(N))
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*sin(x+2*n*pi/N), X(:), 1:N);
% for n = 1:N
%     plot(X(:),Y(:,n), 'linewidth',4);
%     hold all
% end
%
%%% LINE using matrices:
% N = 20;
% set(0,'DefaultAxesColorOrder',tab20c(N))
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*cos(x+2*n*pi/N), X(:), 1:N);
% line(X(:),Y)
%
%% Input and Output Arguments %%
%
%%% Inputs (*=default):
% N = NumericScalar, N>=0, an integer to define the colormap length.
%   = *[], use the length of the current figure's colormap (see COLORMAP).
%
%%% Outputs:
% map = NumericMatrix, size Nx3, a colormap of RGB values between 0 and 1.
%
% See also TAB10 TAB20 TAB20B SET TWILIGHT VIRIDIS LINES COLORMAP PARULA

if nargin<1 || isnumeric(N)&&isequal(N,[])
	N = size(get(gcf,'colormap'),1);
else
	assert(isscalar(N)&&isfinite(N)&&isreal(N),...
		'SC:tab20c:N:NotRealFiniteScalarNumeric',...
		'First argument <N> must be a real finite numeric scalar.')
end
%
hex = ['#3182bd';'#6baed6';'#9ecae1';'#c6dbef';'#e6550d';'#fd8d3c';'#fdae6b';'#fdd0a2';'#31a354';'#74c476';'#a1d99b';'#c7e9c0';'#756bb1';'#9e9ac8';'#bcbddc';'#dadaeb';'#636363';'#969696';'#bdbdbd';'#d9d9d9'];
raw = sscanf(hex.','#%2x%2x%2x',[3,Inf]).';
%
map = raw(1+mod(0:N-1,size(raw,1)),:) / 255;
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tab20c

author: kavli-ntnu


text:
function y = feval(f, x)
%FEVAL   Evaluate a CHEBTECH.
%   Y = FEVAL(F, X) evaluates of the CHEBTECH F at points X.
%
%   If size(F, 2) > 1 then FEVAL returns values in the form [F_1(X), F_2(X),
%   ...], where size(F_k(X)) = size(X).
%
%   Example:
%     f = chebtech2(@(x) 1./( 1 + 25*x.^2 ) );
%     x = linspace(-1, 1, 1000);
%     [xx, yy] = meshgrid(x, x);
%     ff = feval(f, xx + 1i*yy);
%     h = surf(xx, yy, 0*xx, angle(-ff));
%     set(h, 'edgealpha', 0)
%     view(0, 90), shg
%     colormap(hsv)
%
% See also BARY, CLENSHAW, NUDCT.

% Developer note: We use either Clenshaw's algorithm or a nonuniform DCT,
% depending on the dimensions of F and X (determined heuristically).

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

if ( isempty(f) )
    y = [];
    return
end

% Reshape x to be a column vector
[n, m] = size(f);
sizex = size(x);
ndimsx = ndims(x);
x = x(:);

if ( (m > 1) && (ndimsx > 2) )
    error('CHEBFUN:CHEBTECH:feval:evalArrayAtNDArray', ...
        ['Evaluation of a CHEBTECH with more than one column at inputs ' ...
         'with more than two dimensions is not supported.']);
end

if ( n <= 4000 || numel(x) <= 4000 )
    % Evaluate using Clenshaw's algorithm:
    y = f.clenshaw(x, f.coeffs);
else
    % Use fast transform for high degree Chebyshev expansions:
    y = chebfun.ndct(x, f.coeffs);
end

% Reshape the output if possible:
if ( (m == 1) && ( (ndimsx > 2) || (sizex(2) > 1) ) )
    y = reshape(y, sizex);
elseif ( (m > 1) && ( (ndimsx == 2) || (sizex(2) > 1) ) )
    y = reshape(y, sizex(1), m*numel(x)/sizex(1));
end

end


text:
%% polyLineIntersect
% Below is a demonstration of the features of the |polyLineIntersect| function

%%
clear; close all; clc;

%% Syntax
% |[Vn]=polyLineIntersect(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON