NOW/optimizationProblem.m at master · sbaete/NOW · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
classdef optimizationProblem
    %OPTIMIZATIONPROBLEM Store all properties needed to run NOW.
    %   Parameters not specified by the user are default-initialized as follows:
    %
    %   Max gradient = 80 milliTesla/m
    %   Max slew rate = 100 milliTesla/m/milliSecond = 100 T/m/s
    %   Eta (heat dissipation parameter) = 1
    %   Discretization points = 77
    %   Target tensor = eye(3)
    %   Initialguess = 'random'
    %   zeroGradientAtIndex = [], i.e. only at start and end
    %   enforceSymmetry = false;
    %   redoIfFailed = true;
    %   useMaxNorm = false;
    %   doMaxwellComp = true;
    %   MaxwellIndex = 100;
    %   Motion compensation: disabled (when enabled, magnitude unit is s^order / m)


    properties (Access = public)
        targetTensor = eye(3); % Isotropic encoding tensor
        N = 77;
        initialGuess = 'random';
        useMaxNorm = false;
        gMax = 80;
        sMax = 100;
        durationFirstPartRequested = 28;
        durationSecondPartRequested = 22;
        durationZeroGradientRequested = 8;
        eta = 1;
        enforceSymmetry = false;
        redoIfFailed = true;
        name = 'NOW';
        x0 = [];
        doMaxwellComp = true;
        MaxwellIndex = 100;
        MaxFunEval = 1e5;
        MaxIter    = 5e3;
        motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', [])
    end

    properties (SetAccess = private)
        zeroGradientAtIndex = [];
        tolIsotropy = .5e-2;
        tolMaxwell
        signs
        tolSlew
        durationFirstPartActual
        durationZeroGradientActual
        durationSecondPartActual
        totalTimeActual
        dt
        gMaxConstraint
        sMaxConstraint
        integralConstraint
    end

    methods (Access = public)
        function obj = optimizationProblem(varargin)
            if nargin > 0
                settings = varargin{1};

                % Overwrite defaults with user-specified settings
                fieldNames = fieldnames(settings);
                for i = 1:length(fieldNames)
                    eval(['obj.' fieldNames{i} ' = getfield(settings, fieldNames{i});'])
                end
            end

            % Get actual times after discretization
            [obj.durationFirstPartActual, obj.durationZeroGradientActual, obj.durationSecondPartActual, obj.totalTimeActual, obj.zeroGradientAtIndex] = ...
                getActualTimings(obj.durationFirstPartRequested, obj.durationZeroGradientRequested, obj.durationSecondPartRequested, obj.N, obj.enforceSymmetry);


            % Compute private variables
            obj.dt = obj.totalTimeActual/obj.N; %Time step in milliseconds. Division by N instead of N-1 due to half step shift in gradients.
            obj.gMaxConstraint = obj.gMax*obj.dt;
            obj.sMaxConstraint = obj.sMax*obj.dt^2;
            obj.integralConstraint = obj.eta*obj.gMaxConstraint^2*obj.totalTimeActual/obj.dt;

            %% Maxwell compensation
            obj.tolMaxwell = obj.MaxwellIndex/obj.dt; %

            if ~isempty(obj.zeroGradientAtIndex) && obj.doMaxwellComp
                signs = ones(obj.N - 1,1); % Ghost points excluded during opt
                signs(obj.zeroGradientAtIndex(end) + 1:end) = -1;
                obj.signs = signs;

            else
                % Maxwell terms cannot be compensated if no 180 pulses are
                % used. Normally this kind of optimization is intended for
                % a repetition of two identical self-balanced waveforms,
                % but it may be necessary to warn users that single-sided
                % experiments will always incurr some error due to
                % concomitant fields. In practice the "weight" of the
                % optimization with respect to Maxwell terms is removed by
                % setting the sign vector to all zeros.
                obj.doMaxwellComp = false;
                obj.signs = zeros(obj.N - 1,1)+eps; % setting to zero flips out due to sqrt(0)=complex (??)
            end

            %% Motion compensation
            if length(obj.motionCompensation.maxMagnitude) ~= length(obj.motionCompensation.order)
                error('motionCompensation.maxMagnitude must have the same size as motionCompensation.order.')
            end

            if isempty(obj.motionCompensation.maxMagnitude)
                obj.motionCompensation.linear = [];
            else
                % Infer empty motionCompensation.linear from values of
                % motionCompensation.maxMagnitude
                obj.motionCompensation.linear = (obj.motionCompensation.maxMagnitude <= 0);
            end
        end
    end
end