WSINDy_PDE/get_lib_columns.m at master · MathBioCU/WSINDy_PDE · 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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% WSINDy_PDE: compute G and b, listed together in Theta_pdx,
%%%%%%%%%%%% using convolutions over separable test functions via convNDfft
%%%%%%%%%%%%
%%%%%%%%%%%% Copyright 2020, All Rights Reserved
%%%%%%%%%%%% Code by Daniel A. Messenger
%%%%%%%%%%%% For Paper, "Weak SINDy for Partial Differential Equations"
%%%%%%%%%%%% by D. A. Messenger and D. M. Bortz

function Theta_pdx = get_lib_columns(n,lib_list,U_obs,Cfs_x,Cfs_t,m_x,m_t,dx,dt,sub_inds,dim,scales)

Ns = size(U_obs{1});

if isempty(scales)
    scales = ones(n+dim,1);
end

Cfs_ffts = cell(dim,1);
[mm,nn] = size(Cfs_x);
for k=1:dim-1
    Cfs_ffts{k} = [zeros(mm,Ns(k)-nn) (m_x*dx*scales(n+k)).^(-(0:mm-1)').*Cfs_x];
    Cfs_ffts{k} = fft(Cfs_ffts{k},[],2);
end
[mm,nn] = size(Cfs_t);
Cfs_ffts{dim} = [zeros(mm,Ns(dim)-nn) (m_t*dt*scales(n+dim)).^(-(0:mm-1)').*Cfs_t];
Cfs_ffts{dim} = fft(Cfs_ffts{dim},[],2);

Theta_pdx = [];
ind = 1;

while ind<size(lib_list,1)+1
    tags = lib_list(ind,1:n);
    if isreal(tags)
        fcn = (U_obs{1} / scales(1)).^tags(1);
        for k=2:n
            fcn = fcn.*((U_obs{k} / scales(k)).^tags(k));
        end
    else
        ind_freq = find(imag(tags(1:n)));
        freq = sum(imag(tags(1:n)));
        if freq<0
            fcn = sin(abs(freq)*U_obs{ind_freq});
        else
            fcn = cos(abs(freq)*U_obs{ind_freq});
        end
    end

    while all(lib_list(ind,1:n) == tags)
        test_conv_cell = {};
        for k=1:dim-1
            test_conv_cell{k} = Cfs_ffts{k}(lib_list(ind,n+k)+1,:);
        end
        test_conv_cell{dim} = Cfs_ffts{dim}(lib_list(ind,n+dim)+1,:);
        fcn_conv = convNDfft(fcn,test_conv_cell,sub_inds,2);
        Theta_pdx(:,ind) = fcn_conv(:);
        ind = ind+1;
        if ind > size(lib_list,1)
            break
        end
    end
end
end


% function Theta_pdx = get_lib_columns(n,lib_list,U_obs,Cfs_x,Cfs_t,supp_phi_x,supp_phi_t,dx,dt,sub_inds,dim,scales)
%
% Theta_pdx = [];
% ind = 1;
%
% if isempty(scales)
%     while ind<size(lib_list,1)+1
%         tags = lib_list(ind,1:n);
%         if isreal(tags)
%             fcn = U_obs{1}.^tags(1);
%             for k=2:n
%                 fcn = fcn.*(U_obs{k}.^tags(k));
%             end
%         else
%             ind_freq = find(imag(tags(1:n)));
%             freq = sum(imag(tags(1:n)));
%             if freq<0
%                 fcn = sin(abs(freq)*U_obs{ind_freq});
%             else
%                 fcn = cos(abs(freq)*U_obs{ind_freq});
%             end
%         end
%
%         while all(lib_list(ind,1:n) == tags)
%             test_conv_cell = {};
%             for k=1:dim-1
%                 test_conv_cell{k} = Cfs_x(lib_list(ind,n+k)+1,:)' * ((supp_phi_x*dx)^(-lib_list(ind,n+k))) * 1/(2*supp_phi_x+1);
%             end
%             test_conv_cell{dim} = Cfs_t(lib_list(ind,end)+1,:)' * ((supp_phi_t*dt)^(-lib_list(ind,end))*dt) * 1/(2*supp_phi_x+1);
%             fcn_conv = convNDfft(fcn,test_conv_cell,sub_inds,1);
%             Theta_pdx(:,ind) = fcn_conv(:);
%             ind = ind+1;
%             if ind > size(lib_list,1)
%                 break
%             end
%         end
%     end
% else
%     while ind<size(lib_list,1)+1
%         tags = lib_list(ind,1:n);
%         if isreal(tags)
%             fcn = (U_obs{1} / scales(1)).^tags(1);
%             for k=2:n
%                 fcn = fcn.*((U_obs{k} / scales(k)).^tags(k));
%             end
%         else
%             ind_freq = find(imag(tags(1:n)));
%             freq = sum(imag(tags(1:n)));
%             if freq<0
%                 fcn = sin(abs(freq)*U_obs{ind_freq});
%             else
%                 fcn = cos(abs(freq)*U_obs{ind_freq});
%             end
%         end
%
%         while all(lib_list(ind,1:n) == tags)
%             test_conv_cell = {};
%             for k=1:dim-1
%                 test_conv_cell{k} = Cfs_x(lib_list(ind,n+k)+1,:)' * ((supp_phi_x*scales(n+k)*dx)^(-lib_list(ind,n+k)) * 1/(2*supp_phi_x+1));
%             end
%             test_conv_cell{dim} = Cfs_t(lib_list(ind,end)+1,:)' * ((supp_phi_t*scales(end)*dt)^(-lib_list(ind,end)) * 1/(2*supp_phi_t+1));
%             fcn_conv = convNDfft(fcn,test_conv_cell,sub_inds,1);
%             Theta_pdx(:,ind) = fcn_conv(:);
%             ind = ind+1;
%             if ind > size(lib_list,1)
%                 break
%             end
%         end
%     end
% end