forked from ruthamey/slipBERI
-
Notifications
You must be signed in to change notification settings - Fork 0
/
slipBERI.m
2817 lines (2368 loc) · 175 KB
/
slipBERI.m
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function [ ] = slipBERI( fault, data, invert, priors, elastic_params, display, housekeeping )
%
% slipBERI is a code to invert for distribued earthquake slip,
% incorporating fractal properties. It does so using Bayesian methods.
%
% Several lines of evidence suggest that earthquake slip should show
% fractal properties. This includes fault surface roughness, measurements
% of surface slip during an earthquake and analysis of published slip
% solutions. In light of this fractal properties are incorporated into the
% slip inversion as a prior, through the von Karman autocorrelation. Markov
% chain Monte Carlo sampling is used to sample parameter space, using the
% Metropolis-Hastings rejection.
%
% The inputs to this function are seven structures, containing all the
% details required to run, display and save the results of the inversion.
% The provided code 'make_structures_required_for_slipBERI' should be
% edited as appropriate and then run to set-up these initial structures.
% Or the variables can be changed from the terminal:
% structure.name_of_variable = <new_value>
%
% All of the details of the inputs are listed at the start of the slipBERI
% code, or for a more reader-friendly version please see the google docs
% help here:
% https://docs.google.com/document/d/1cUXLRxN-oB8Q8kGOueq2c-Zxr3W1vDgWpGUw1MAEx5s/edit?usp=sharing
%
% This function relies on many scripts, found in 'slipBERImbin'. Each script should
% give full details of what it requires, does, and outputs, found by typing
% 'help <name_of_script>'
%
% To run this code, simply enter into the matlab terminal:
% slipBERI( fault, data, invert, priors, elastic_params, display, housekeeping )
%
% Written by Ruth Amey (rmja) 18th-Nov-2014 onwards...
% With Andy Hooper & some scripts from GJ Funning, TJ Wright, D Bekaert, T Ingleby
%
% When using this code please cite:
% Amey, R. M. J., Hooper, A., & Walters, R. J. (2018). A Bayesian method
% for incorporating self‐similarity into earthquake slip inversions.
% Journal of Geophysical Research: Solid Earth, 123, 6052–6071.
% or
% Amey, R.M.J., Hooper, A. and Morishita, Y. Going to Any Lengths: Solving
% for Fault Size and Fractal Slip for the 2016, Mw 6.2 Central Tottori
% Earthquake, Japan, using a Trans-dimensional Inversion Scheme, Journal of
% Geophysical Research: Solid Earth
% *********** % ************ % ********** % ********** % ********** % *****
% The inputs are seven structures containing this information:
%
% INPUTS:
% 'fault' is a structure relating to the description of the fault upon which you'd like to invert:
% fault_descriptor_file = String. Name of text file of: [strike, dip, rake, up dip projection of CENTRE of fault plane, up dip projection of CENTRE of fault plane, length (km), top depth (km), bottom depth (km), patches along strike, patches down dip, smoothingidentifyer]. where fault strands with the same smoothingidentifyer will be smoothed across, if you select to smooth across fault strands at all. If you please.
% fault_coordinate_unit = String. 'utm' or 'long/lat' to indicate whether the coordinates used in fault.fault_desriptor_file are in given in UTM or long/lat.
%
% 'data' is a structure containing details of the data that you're inverting:
% InSAR_datafile = Text file of [x (long), y (lat), observed LOS displacement, (E, N, up) component of LOS vector]. Make sure you've downsampled and added pointing vector. LOS vector convention is positive looking at the satellite. If you have no InSAR data this MUST be 'none'. Locations translated to UTMx or a local latitude, if spanning multiple UTM zones.
% InSAR_coordinate_unit = 'utm' or 'long/lat'. Everything will be converted to UTM as necessary.
% varcovar_details = String. this is a text file with the [sill, variogram_range] from the covariogram from an undeforming region. Can be calculated in advance using andy's code variogram.m; In older versions this was: Name of a text file containing the outputsfrom cvdcalc: [maxium covariance in metres^2, alpha from expcos function, beta from expcos function, eb_r from ebessel function, eb_w from ebessel function, either 'expcos' or 'ebessel' depending on which is the best function);].
% quadtree_n_points = String. Name of a text file with the number of pixels combined into each datapoint, if using quadtree. This is used to calculate the var-covar matrix. If not using quadtree this is a text file with '1's in each row, the same length as the InSAR datafile.
% GPS_datafile_2d = GPS file of 2d surface displacement. Set upto be the same as InSAR. Text file of 7 columns [x (long), y (lat), observed LOS displacement (m), LOS vector E, LOS vecor N, LOS vector up, sigma] where sigma is the standard deviation (mm for GPS, m for atolls). Each GPS station will have two rows: one for E LOS vector, one for N LOS vector. The sigma must be for E, N or up appropriately depending on which value you have. If you have no GPS data this MUST be 'none'. Locations translated to UTMx or a local latitude, if spanning multiple UTM zones.
% GPS_datafile_3d = GPS file of 3d surface displacement. Set upto be the same as InSAR. Text file of 7 columns [x (long), y (lat), observed LOS displacement (m), LOS vector E, LOS vecor N, LOS vector up, sigma] where sigma is the standard deviation (mm for GPS, m for atolls). Each GPS station will have three rows: one for E LOS vector, one for N LOS vector, one for up LOS vector. The sigma must be for E, N or up appropriately depending on which value you have. If you have no GPS data this MUST be 'none'. Locations translated to UTMx or a local latitude, if spanning multiple UTM zones.
% GPS_coordinate_unit = String. 'lat/long' or 'UTM'. Everything will be converted to UTM as necessary.
% atolls_datafile = String. Name text file with uplift only data in.
% atolls_coordinate_unit = String. 'utm' or 'long/lat'
% weight_InSAR = Number. Relative weighting of InSAR data. If not weighting datasets then this must be 1. If this is a number other than 1 then the InSAR var-covar matrix will be divided by this value.
% weight_GPS = Number. Relative weighting of GPS data. If not weighting datasets then this must be 1. If this is a number other than 1 then the GPS var-covar matrix will be divided by this value.
% weight_atolls = Number. Relative weighting of atoll dataset.
% UTMzone = Number (ignoring letter). The UTM zone e.g. 10
% origin = A long/lat in the centre of your area [-124 38]. This is only used if your data spans multiple UTM zones and in that case everything is converted to a local coordinate system, using the origin.
% EQ_epicenter = Number. Matrix of [lat, long, depth]. Depth in km.
% seismic_moment = Number. Only used if invert.regularise_moment == 'yes'. Estimate of seismic moment e.g. from USGS in Nm
% moment_std = Number. Only used if invert.regularise_moment == 'yes'. Estimate of seismic moment e.g. variation of USGS solutions
%
% 'invert' is a structure with details of how you wish the inversion to be performed:
% quickcheck = String. 'yes' or 'no'. Before commencing the inversion, your data will be plotted for you to click on to say that it's okay.
% inversion_type = String. 'least_squares' (shame on you), or 'bayesian'
% iterations = Number. Number of iterations, only relevant for Bayesian inversions.
% smoothing = String. 'None' gives no smoothing, 'laplacian' does Laplacian smoothing, 'VK' does von Karman smoothing.
% smooth_across_fault_strands = String. 'yes' smoothes across different fault strands as if it were one fault. 'no' treats each fault as a different earthquakes. [] means I'll calculate which option you should choose based on the fault geometry .
% slip_initial = Number. Estimate of slip values - all slip patches will be assigned this amount of slip initially, with a bit of noise added.
% step_size = Value. Maximum stepsize, plus or minus, allowed in your random walk when generating a slip_trial, for Bayesian inversions. This is used initially and then is adjusted through the iteration llh2local. For 'invert.slip_prior = boxcar', step_size is in metres. For 'invert.slip_perior = logarithmic' it's like a percentage of slip, or something.
% variable_rake = String. 'yes' or 'no'
% solve_for_InSAR_offset: String. 'no' or 'yes'
% solve_for_InSAR_ramp: String. 'no' or 'yes'
% regularise_moment = String. 'yes' or 'no'. This adds an M0 prior likelihood - a normal distribution, with mean and std given in 'data' structure.
% alpha2_initial = Number. One number for each fault strand, of initial alpha^2 (variance) value. This is adjusted through the iteration.
% alpha2_step_size = Number. Initial step size of alpha2_modelparameter, which is then adjusted through the inversion.
% probability_target_initial = Number. Initial probability target - if we add half the step size to a parameter, this is the probability we're aiming the perturbation to make. This is adjusted through the iteration process: if the rejection rate is too high, then the probability decreases to try to decrease step sizes. % JUST FOR NOW while I'm sorting out sensitivity
% solve_for_beta = String. 'yes' or 'no' - beta is a hyperparameter on the data.
% beta_initial = Number. Starting value of hyperparameter that acts on the var-covariance matrix. If this is set to 1 and beta_step_size is set to 0 then beta will stay at 1 throughout the whole inversion.
% beta_step_size = Number. Starting value of hyperparmaeter stepsize. Default is 0.
% simulated_annealing_start = String. 'yes' or 'no'. slipBERI will perform an initial simulated annealing inversion, to use as the starting parameter for the Bayesian inversion.
% solve_for_fault_size: String. 'yes' or 'no'. NOTE: THIS ONLY WORKS WITH ONE FAULT STRAND, currently.
% add_correlation_matrix_stabiliser = String. 'yes' or 'no'. Add a small term (akin to a nugget) to the diagonal of the von Karman matrix sigma_s. Unlikely to be necessary unless solving for the fault size.
% load_old_MCMC_chain = String. Either a name of a saved file, or 'no'. If this is a name of a file, then the MCMC chain will continue from the past max-likelidhood solution, with same step sizes, alpha^2, and probability target. BE CAREFUL using this if you're loading a chain with DIFFERENT inversion parameters. some things may not work.
%
% 'priors' is a structure containing the priors
% slip_prior = String. 'boxcar' or 'gaussian' or 'logarithmic', for Bayesian inversions. Note that if you use logarithmic, you will probably need to increase your number of iterations, since there are more rejections, and also maybe decrease your step size.
% min_slip = Number. Minimum value of slip allowed, in metres.
% max_slip = Number. Maximum value of slip allowed, in metres.
% predominant_faulting_style: String in curly brackets. Either {'ss'} or {'ds'}. This is used in von Karman smoothing to calculate the along-strike and down-dip correlation lengths.
% min_rake = Number. Minimum value of permitted rake. Rake 0 = left-lateral strike-slip. Rake 180 = right-lateral strike-slip. Rake 90 = thrust. Rake -90 (or 270) = normal. Can be a matrix of one value per each fault_strand_for_smoothing, or if just one value used then same value is used across all fault strands.
% max_rake = Number. Maximum value of permitted rake. Rake 0 = left-lateral strike-slip. Rake 180 = right-lateral strike-slip. Rake 90 = thrust. Rake -90 (or 270) = normal. Can be a matrix of one value per each fault_strand_for_smoothing, or if just one value used then same value is used across all fault strands.
% alpha2_prior = String. 'logarithmic' or [], depending on if you want to use a logarithmic prior on alpha2 or not.
% max_alpha2 = Number. Maximum permitted alpha^2 hyperparameter. Note that if using logarithmic prior this is the maximum alpha^2 value, not maximum 10^(alpha^2) value.
% min_alpha2 = Number. Minimum permitted alpha^2 hyperparameter. Note that if using logarithmic prior this is the minimum alpha^2 value, not maximum 10^(alpha^2) value.
% alpha2_flag = String. Either 'bothsame' if you wish to use the same alpha2 hyperparameter on multiple fault strands or [] if not.
% min_dip = Number. Minimum value of dip permitted. slipBERI CANNOT CURRENTLY SOLVE FOR DIP -- to be added soon. If there are multiple fault strands, can enter one value per fault strand. If only one number is entered this is used for all fault strands.
% max_dip = Number. Maximum value of dip permitted. slipBERI CANNOT CURRENTLY SOLVE FOR DIP -- to be added soon. If there are multiple fault strands, can enter one value per fault strand. If only one number is entered this is used for all fault strands.
% max_offset = Number. Maximum permitted offset if invert.solve_forInSAR_offset = 'yes'.
% min_offset = Number. Minimum permitted offset if invert.solve_forInSAR_offset = 'yes'.
% min_beta = Number. Minimum value of permitted beta. The default is to not solve for beta, so this is only used if invert.solve_for_beta = 'yes'.
% max_beta = Number. Maximum value of permitted beta. The default is to not solve for beta, so this is only used if invert.solve_for_beta = 'yes'.
% min_circharm_coeffs = Number. Minimum permitted circharm coefficient. This is related to the circular harmonics terms if solving for fault size, so this is only used if invert.solve_for_fault_size = 'yes'.
% max_circharm_coeffs = Number. Maximum permitted circharm coefficient. This is related to the circular harmonics terms if solving for fault size, so this is only used if invert.solve_for_fault_size = 'yes'.
% min_circharm_phi = Number. Minimum permitted circharm rotation (radians). This is related to the circular harmonics terms if solving for fault size, so this is only used if invert.solve_for_fault_size = 'yes'.
% max_circharm_phi = Number. Maximum permitted circharm rotation (radians). This is related to the circular harmonics terms if solving for fault size, so this is only used if invert.solve_for_fault_size = 'yes'.
% min_circharm_center = Number. Maximum permitted circharm center location (m). This is related to the circular harmonics terms if solving for fault size, so this is only used if invert.solve_for_fault_size = 'yes'.
% max_circharm_center = Number. Maximum permitted circharm center location (m). This is related to the circular harmonics terms if solving for fault size, so this is only used if invert.solve_for_fault_size = 'yes'.
%
% 'elastic_params' is a structure containing the Lame elastic constants whic are required for kernel calculation (Okada 1985)
% lambda = 3.23e10. Number. Value of Lame's first parameter.
% mu_okada = 3.23e10. Number. Value of Lame's second parameter, shear modulus.
%
% 'display' is a structure which details how you would like your data to be displayed.
% plotmean = String, 'yes' or 'no'. whether you'd like to plot the mean and standard deviation or not.
% plotmode = String, 'yes' or 'no'. whether you'd like to plot the mode and standard deviation or not.
% plotmedian = String, 'yes' or 'no'. whether you'd like to plot the median and standard deviation or not.
% plotallsips = String, 'yes' or 'no'. whether you'd like one plot with the mean, mode, median, max likelihood (and true, if if in testing mode).
% plotprob = String, 'yes' or 'no'. whether you'd like to plot how the probability changes with the number of iterations.
% plothists = String. how many histograms you'd like to display. 'plothistall' means plot histograms for all the slip patches, 'plothistsample' means select some randomly to display
% plotsurfacedisp = String, 'yes' or 'no'. Display surface displacement - arrows for GPS, colour for InSAR LOS.
% plotmarginalPDFs = String, 'yes' or 'no'. this chooses the six patches with the highest slip and plots the marginal PDFs for this.
% plot_resolution_matrix = String. 'yes' or 'no'.
% plot3d = String. 'yes' or 'no' for whether to plot the fault geometry with slip mode.
% plotMAP = String. 'yes' or 'no' for whether to plot the MAP (max a posteriori) and standard deviation or not
% calc_confidence = String. 'yes' or 'no' for whether you want to plot the 95% confidence intervals for each patch.
%
% 'housekeeping' is important for keeping your files (as well as your mind, body and soul) tidy:
% save_name = String. name you'd like to call your data run when you save. Note that slipBERI automatically appends it with your inversion style
%
% *********** % ************ % ********** % ********** % ********** % *****
%% CHANGE HISTORY - for the backseat drivers amongst you ******************
% 'Date, initial, details' if you please. Let's keep it tidy in here.
% deets: rmja ([email protected])
% 18-nov-2014 rmja Started writing
% 26-feb-2015 rmja Code actually ran the whole way through
% 08-may-2015 rmja Moved from a script to a function
% 29-jul-2015 rmja Fixed plotting InSAR, GPS and residuals
% 31-jul-2015 rmja Added variable rake
% 18-aug-2015 rmja Added converting GPS from lat/long to UTM
% 24-aug-2015 rmja changed 'logprob_VK_temp' to 'logprior_VK_temp'; 'logprob_VK_trial' to 'logprior_VK_trial'; and 'logprob_VK_curr' to 'log_prior_VK_curr'.
% 7-sep-2015 rmja Solving for alpha_squared as a hyperparameter
% 28-sep-2015 rmja Fixed sensitivity - aiming for an ideal probability perturbation, which is also solved for
% 7-oct-2015 rmja Add NaN checks
% 8-oct-2015 rmja Added zero padding option on left, right, bottom edges of fault.
% 22-oct-2015 rmja Corrected calculating RMS to WRMS
% 27-oct-2015 rmja Separated 2D and 3D GPS
% 2-nov-2015 rmja Can now solve for multiple fault strands
% 9-nov-2015 rmja Remove points a mindist from the fault
% 5-jan-2016 rmja Added Laplacian smoothing into MCMC. changed logprior_VK_curr to logprior_curr, logprior_VK_temp to logprior_temp and logprior_VK_trial to logprior_trial
% 5-jan-2016 rmja Can choose whether to smooth across multiple fault strands or not
% 14-jan-2016 rmja Added in M0 regularisation
% 10-feb-2016 rmja Solve for along-strike and down-dip Vk correlation lengths
% 26-feb-2016 rmja Solve for dip
% 16-jun-2016 rmja Hacked process_faultdata_center a LOT so that when I plot results using imagesc it always plots the northerly and westerly patch in the top left.
% 4-jul-2016 rmja Correcting weighting different datasets
% 8-jul-2016 rmja Sorted out utm/latlong/local coordinates
% 13-jul-2016 rmja Added in atolls as a dataset
% 14-jul-2016 rmja Added in solving for InSAR offset
% 4-aug-2016 rmja Fix solving for multiple faults under the surface (ramp flat ramp)
% 9-dec-2016 rmja Added in solving for fault length
% 1-feb-2017 rmja Changed varcovar calculation to calculating variogram
% 13-feb-2017 rmja Fixed smoothing so you can identify which fault strands you want to smooth across and which you don't (whereas previously was smooth across all or smooth across none)
% 14-feb-2017 rmja Fixed solving for patches being on/off so you can now use it with multiple fault strands
% 4-aug-2017 rmja So you can choose different min/max rakes on different fault strands
% 23-aug-2017 rmja Solve for fault size with circular harmonics
% 25-aug-2017 rmja Ensemble approach to taking step sizes (Goodman and Weare, 2010)
% 15-mar-2018 rmja Parallelising ensemble move (Foreman-Mackey et al., 2013).
% 27-mar-2018 rmja Tidied code and fixed all the archaeic nonsense, lurking in subfunctions
% 22-may-2018 rmja Solving for correlation lenghts a_as and a_dd from their distributions, when solving for fault size.
% 15-jan-2019 rmja Tidied significantly for github
% 20-jan-2019 rmja Added solving for ramp
%% Sort out some matrices.
% These are half-coded options that I don't want to advertise and make
% available to the general user. However if you've found these and are
% interested, then these are options that you could debug and test. Most of
% them are set up for one fault strand, it just might be inefficient or
% doesn't work for every option.
% Some inverting parameters that I did not finish
invert.solve_for_dip = 'no'; % String. 'yes' or 'no' to solve for dip. Not debugged, very slow, may or may not work.
invert.select_patches = 'no'; % String. The name of the matlab file that contains the numbers of which slip patches you wish to be 'on' during this inversion.
invert.ensemble_sampling = 'no'; % String. 'yes' or 'no' for whether you want to explore parameter space using ensemble sampling, rather than adding a random number to the model parameters.
invert.ensemble_move_type = []; % String. If using ensemble sampling you can choose whether you wish to use the 'stretch' or 'walk' method of doing so.
invert.solve_for_correlation_length = 'no'; % String. 'yes' or 'no'. This is the option to solve for the along-strike and down-dip correlation terms for VK smoothed solution
% Testing model, used to make the synthetic tests in Amey et al. 2018.
% They require a very specific set-up, which I have not provided, so I have
% not made it available to the user
testing.testing_mode = 'no'; % String. 'yes' or 'no'. Note that if set to 'yes' any file in data.InSAR_datafile or data.GPS_datafile will be ignored and instead data in testing.making_model will be used.
testing.making_model = []; % String. File of saved data from 'making_model.m'
testing.var = []; % Number. Imaginary values of variance for dataset. Assume no covariance.
testing.add_noise = []; % Number. 0 = no noise, 1 = 1sigma noise, 2 = 2sigma-noise, etc etc
%% CHECK PATHS
% Now on the help file I get the user to add their path, rather than
% putting it here
%addpath(genpath('/nfs/see-fs-01_users/eermja/Documents/MATLAB/slipBERImbin'));
profile on
%% Let's check our inputs are all fine-and-dandy before we proceed
run('error_messages.m');
%% FUN WITH FUNCTIONS - sit back, relax, let the magic begin **************
tic
% *********** % ************ % ********** % ********** % ********** % *****
% Introductions
% *********** % ************ % ********** % ********** % ********** % *****
disp(' ');
disp('Now running slipBERI - the distributed slip inversion code with fractal regularisation');
disp('by Ruth Amey and Andy Hooper (COMET, University of Leeds)');
disp('Check for updates at https://github.com/ruthamey/slipBERI');
disp(' ')
rng('shuffle') % Fun fact - if you don't set your rng to shuffle then every time you open matlab fresh it'll give you the same random numbers in a row. Try it - I bet you £50 if you open a new matlab and type 'rand' you'll get 0.8147
% *********** % ************ % ********** % ********** % ********** % *****
% Look at the data
% *********** % ************ % ********** % ********** % ********** % *****
use_local_coordinate_system = 'no'; % by default use utmx utmy, not a local coordinate system
if strcmp(testing.testing_mode, 'yes') == 1
load(testing.making_model, 'locs');
n_datasets = 1;
los_vector = ones(1,length(locs));
data.GPS_datafile = 'none';
first_InSAR_scene_numbers = [];
InSAR_identifyer = [];
n_InSAR_scenes = [];
else
n_datasets = 0;
% Sort out InSAR - data, locations, line of sight vector
if strcmp(data.InSAR_datafile, 'none') ~= 1 % if it's TRUE that there is NOT 'none' InSAR data. i.e. if there IS an insar datafile
n_datasets = n_datasets + 1;
d_InSAR = [];
locs_InSAR = [];
los_vector_InSAR = [];
counter = 1;
InSAR_identifyer = [];
n_InSAR_scenes = size(data.InSAR_datafile,2);
first_InSAR_scene_numbers = 1;
for p = 1:n_InSAR_scenes
[ locs_tmp, d_tmp, los_vector_tmp ] = read_InSAR( cell2mat(data.InSAR_datafile(p)));
d_InSAR = [d_InSAR, d_tmp];
locs_InSAR = [locs_InSAR, locs_tmp];
los_vector_InSAR = [los_vector_InSAR, los_vector_tmp];
InSAR_identifyer = [InSAR_identifyer; counter*ones(length(d_tmp),1)];
first_InSAR_scene_numbers = [first_InSAR_scene_numbers; (length(d_tmp)+first_InSAR_scene_numbers(end))];
counter=counter+1;
locs_tmp=[];
los_vector_tmp=[];
d_tmp=[];
end
first_InSAR_scene_numbers(end) = [];
if strcmp(data.InSAR_coordinate_unit, 'long/lat') == 1
locs_InSAR_latlong = locs_InSAR;
[locs_InSAR(1,:),locs_InSAR(2,:)]=ll2utm(locs_InSAR(2,:),locs_InSAR(1,:));
elseif strcmp(data.InSAR_coordinate_unit,'utm') == 1
locs_InSAR = locs_InSAR*1000;
[locs_InSAR_latlong(2,:), locs_InSAR_latlong(1,:)] = utm2ll( locs_InSAR(1,:), locs_InSAR(2,:), repmat(data.UTMzone, 1, length(locs_InSAR)));
end
if locs_InSAR_latlong(1,:) < 0 % want long from
locs_InSAR_latlong(1,:) = locs_InSAR_latlong(1,:) + 360;
end
clear d_tmp
else
d_InSAR = [];
locs_InSAR = [];
n_InSAR_scenes = [];
first_InSAR_scene_numbers = [];
InSAR_identifyer = [];
end
% Sort out GPS - data, locations, line of sight vector
if strcmp(data.GPS_datafile_2d, 'none') + strcmp(data.GPS_datafile_3d, 'none') ~= 2 % if it's TRUE that there is NOT 'none' GPS data. i.e. if there IS an GPS datafile
n_datasets = n_datasets + 1; % GPS just counts as one dataset, even if 2d and 3d. yup.
if strcmp(data.GPS_datafile_2d, 'none') ~= 1 % 2D GPS
[ locs_GPS_2d, d_GPS_2d, los_vector_GPS_2d, sigma_gps_2d ] = read_GPS( data.GPS_datafile_2d );
d_GPS_e_2d = d_GPS_2d(1:2:length(d_GPS_2d));
d_GPS_n_2d = d_GPS_2d(2:2:length(d_GPS_2d));
d_GPS_up_2d = [];
else
d_GPS_2d = [];
d_GPS_e_2d = [];
d_GPS_n_2d = [];
d_GPS_up_2d = [];
locs_GPS_2d = [];
los_vector_GPS_2d = [];
sigma_gps_2d = [];
end
if strcmp(data.GPS_datafile_3d, 'none') ~=1 % 3D GPS
[ locs_GPS_3d, d_GPS_3d, los_vector_GPS_3d, sigma_gps_3d ] = read_GPS( data.GPS_datafile_3d );
d_GPS_e_3d = d_GPS_3d(1:3:length(d_GPS_3d));
d_GPS_n_3d = d_GPS_3d(2:3:length(d_GPS_3d));
d_GPS_up_3d = d_GPS_3d(3:3:length(d_GPS_3d));
else
d_GPS_3d = [];
d_GPS_e_3d = [];
d_GPS_n_3d = [];
d_GPS_up_3d = [];
locs_GPS_3d = [];
los_vector_GPS_3d = [];
sigma_gps_3d = [];
end
d_GPS = [ d_GPS_2d, d_GPS_3d];
d_GPS_e = [ d_GPS_e_2d, d_GPS_e_3d];
d_GPS_n = [ d_GPS_n_2d, d_GPS_n_3d];
d_GPS_up = [ d_GPS_up_2d, d_GPS_up_3d];
locs_GPS = [ locs_GPS_2d,locs_GPS_3d ];
los_vector_GPS = [los_vector_GPS_2d ,los_vector_GPS_3d ];
sigma_gps = [ sigma_gps_2d ; sigma_gps_3d ];
data.GPS_datafile = 'yes';
if strcmp(data.GPS_coordinate_unit, 'long/lat') == 1
[locs_GPS(1,:),locs_GPS(2,:), ~ ] = ll2utm(locs_GPS(2,:), locs_GPS(1,:));
if strcmp(data.GPS_datafile_2d, 'none') ~= 1
[locs_GPS_2d(1,:),locs_GPS_2d(2,:), ~ ] = ll2utm(locs_GPS_2d(2,:), locs_GPS_2d(1,:));
end
if strcmp(data.GPS_datafile_3d, 'none') ~=1
[locs_GPS_3d(1,:),locs_GPS_3d(2,:), ~ ] = ll2utm(locs_GPS_3d(2,:), locs_GPS_3d(1,:));
end
else
%locs_GPS = locs_GPS*1000;
end
locs_GPS_unique(1,:) = unique(locs_GPS(1,:), 'stable');
locs_GPS_unique(2,:) = unique(locs_GPS(2,:), 'stable');
else
data.GPS_datafile = 'none';
d_GPS = [];
locs_GPS = [];
los_vector_GPS = [];
end
% Sort out atolls, if you have atolls
if strcmp(data.atolls_datafile, 'none') ~= 1
[ locs_atolls, d_atolls, los_vector_atolls, sigma_atolls ] = read_GPS( data.atolls_datafile );
% Convert from latlong and name something different
if strcmp(data.atolls_coordinate_unit,'long/lat') == 1
locs_atolls_latlong = locs_atolls;
[locs_atolls(1,:),locs_atolls(2,:), ~ ] = ll2utm(locs_atolls(2,:), locs_atolls(1,:));
end
else
d_atolls = [];
locs_atolls = [];
los_vector_atolls = [];
end
% % see if you need to convert to a local coordinate system or keep as utm
%
% if strcmp(data.GPS_datafile, 'none') ~= 1 && strcmp(data.GPS_coordinate_unit, 'long/lat') ==1
% [~,~, zones_GPS ] = ll2utm(locs_GPS(2,:), locs_GPS(1,:));
% else
% zones_GPS = [];
% end
%
% if strcmp(data.InSAR_datafile, 'none') ~= 1
% [~,~, zones_InSAR ] = ll2utm(locs_InSAR_latlong(2,:), locs_InSAR_latlong(1,:));
% else
% zones_InSAR = [];
% end
%
% if strcmp(data.atolls_datafile, 'none') ~= 1
% if strcmp(data.atolls_coordinate_unit, 'long/lat') == 1
% [~,~, zones_atolls ] = ll2utm(locs_atolls_latlong(2,:), locs_atolls_latlong(1,:));
% else
% zones_atolls = [];
% end
% else
% zones_atolls = [];
% end
%
% zones_test = [zones_GPS, zones_InSAR, zones_atolls];
%
% if any(diff(zones_test) ~= 0 )
% disp('having commented out converting to local coordinate system if UTM zones are different')
% use_local_coordinate_system = 'yes';
% disp('Converting GPS and InSAR to local coordinate sytem...')
%
% if strcmp(data.GPS_datafile, 'none') ~= 1
% llh_GPS = [locs_GPS(1,:); locs_GPS(2,:); zeros(1,length(locs_GPS))];
% xy_GPS = llh2local(llh_GPS, data.origin);
% locs_GPS(1,:) = xy_GPS(1,:)*1000; % convert to meters
% locs_GPS(2,:) = xy_GPS(2,:)*1000; % convert to meters
% end
%
% if strcmp(data.InSAR_datafile, 'none') ~= 1
% llh_InSAR = [locs_InSAR(1,:); locs_InSAR(2,:); zeros(1,length(locs_InSAR(1,:)))];
% xy_InSAR = llh2local(llh_InSAR, data.origin);
% locs_InSAR(1,:) = xy_InSAR(1,:)*1000; % convert to meters
% locs_InSAR(2,:) = xy_InSAR(2,:)*1000; % convert to meters
% end
%
% end
% Combine all data, locations, line of sights
d = [d_InSAR, d_GPS, d_atolls];
locs = [locs_InSAR, locs_GPS, locs_atolls];
los_vector = [los_vector_InSAR, los_vector_GPS, los_vector_atolls];
d = d';
n_data = length(d);
d_identifyer = [ 1* ones(length(d_InSAR),1); 2*ones(length(d_GPS),1); 3*ones(length(d_atolls),1)];
end
% *********** % ************ % ********** % ********** % ********** % *****
% Sort out the fault plane
% *********** % ************ % ********** % ********** % ********** % *****
% Create dislocation model for each patch
[disloc_model, spatial_model1, spatial_model2, spatial_model3, spatial_model4, fault_coords, ~, n_along_strike_patches, n_down_dip_patches, fault_length, fault_width,n_fault_strands, strike, fault_strand_togetherness ] =...
process_faultdata_centre(fault, invert, use_local_coordinate_system, data.origin, testing); % if you're padding the edges with zeros, THESE EXTRA PATCHES GET ADDED HERE
spatial_model2column = reshape(spatial_model2, [],1);
spatial_model3column = reshape(spatial_model3, [],1);
% Since we're not solving for dip at the moment...
priors.min_dip = disloc_model(4,1);
priors.max_dip = disloc_model(4,1);
% Since we're not solving for dip at the moment...
% %Save for GMT
fault_coords_latlong = [];
for i = 1:n_fault_strands
[LAT,LON]=utm2ll([fault_coords(i,1), fault_coords(i,3)],[fault_coords(i,2), fault_coords(i,4)],data.UTMzone);
fault_coords_latlong = [fault_coords_latlong; (LON)', LAT'];
end
dlmwrite('fault_coords_latlong.gmt', fault_coords_latlong);
% Start with some easy thing
%total_n_slip_patches = n_along_strike_patches' * n_down_dip_patches; % this works even if you have multiple fault strands
total_n_slip_patches = size(disloc_model, 2);
n_slip_patches_on_each_fault_strand = n_along_strike_patches .* n_down_dip_patches; % one value per fault strand
first_patch_in_strand(1) = 0;
first_patch_in_strand(2:(n_fault_strands+1),1) = cumsum(n_slip_patches_on_each_fault_strand);
first_patch_in_strand = first_patch_in_strand+1;
first_patch_in_strand(n_fault_strands+1) = []; % so that you know which number slip patch is the first one in each new strand. comes in handy
last_patch_in_strand = cumsum(n_slip_patches_on_each_fault_strand); % so that you know which number slip patch is the last one in each new strand. comes in handy
onoffidentifyer = ones(total_n_slip_patches,1);
along_strike_sep_dist = zeros(n_fault_strands,1);
down_dip_sep_dist = zeros(n_fault_strands,1);
fault_strand_identifyer = [];
fault_strand_identifyer_for_smoothing = [];
difference_between_faultssegments_and_smoothingfaultsegments = diff(fault_strand_togetherness);
difference_between_faultssegments_and_smoothingfaultsegments = [ 0;difference_between_faultssegments_and_smoothingfaultsegments];
counter = 1;
for j = 1 : n_fault_strands
along_strike_sep_dist(j,1) = spatial_model2(first_patch_in_strand(j)); % one number for each strand. ASSUMING ALL FAULT PATCHES ARE THE SAME SIZE IN ONE FAULT STRAND. otherwise I can make a matrix of separation distances and separate terms for each fault patch
down_dip_sep_dist(j,1) = spatial_model3(first_patch_in_strand(j)); % one number for each strand. ASSUMING ALL FAULT PATCHES ARE THE SAME SIZE IN ONE FAULT STRAND.
identifyer = j * ones( n_slip_patches_on_each_fault_strand(j),1);
secondidentifyer = counter * ones( n_slip_patches_on_each_fault_strand(j),1);
fault_strand_identifyer = [fault_strand_identifyer; identifyer ];
if difference_between_faultssegments_and_smoothingfaultsegments(j)==1 % then this is considered a new strand to the last strand
counter = counter + 1;
secondidentifyer = counter * ones( n_slip_patches_on_each_fault_strand(j),1);
fault_strand_identifyer_for_smoothing = [fault_strand_identifyer_for_smoothing; secondidentifyer ]; % one value for each fault patch
elseif difference_between_faultssegments_and_smoothingfaultsegments(j)==0 % then this is considered the same strand as the last strand
fault_strand_identifyer_for_smoothing = [fault_strand_identifyer_for_smoothing; secondidentifyer ]; % one value for each fault patch
end
end
clear identifyer
if any(disloc_model(8,first_patch_in_strand)~=0) % if any of the faults start below ground then we need to sum all the down dip patches, and they'll all have the same number of along strike patches
total_n_along_strike_patches = n_along_strike_patches(1); % this is necessary because sometimes we have multiple strands laterally (e.g. napa) and sometimes we have them going down depth (e.g. nepal)
total_n_down_dip_patches = sum(n_down_dip_patches);
else % if all the faults start at the surface but there are a few extending laterally then we need to sum all the alongstrike patches, and they'll all have the same number of downdip patches
total_n_along_strike_patches = sum(n_along_strike_patches);
total_n_down_dip_patches = n_down_dip_patches(1);
end
% Multiple fault strands...............................................
if strcmp(invert.smooth_across_fault_strands, []) ==1 % if they left it to me to decide whether to treat fault strands differently or not.
if any(abs(diff(strike))) > 18 % austin elliott et al 2015 - bends greater than 18 degrees terminate rupture, so should treat as different faults
invert.smooth_across_fault_strands = 'yes';
else
invert.smooth_across_fault_strands = 'no';
end
end
if strcmp(invert.smooth_across_fault_strands, 'yes') == 1 % note that if you do this then only the fault strands with the SAME IDENTIFYER IN the fault geometry file will be smoothed across
%n_fault_strands_for_smoothing = 1;
n_fault_strands_for_smoothing = length(unique(fault_strand_togetherness));
%fault_strand_identifyer_for_smoothing = []; % we made this assuming all the different fault strands would be separate. so we'll ignore that now
difference_between_faultssegments_and_smoothingfaultsegments = diff(fault_strand_togetherness);
for r = 1: n_fault_strands_for_smoothing % this works out how many fault strands you wanna have
fault_length_for_smoothing(r,1) = sum(fault_length(fault_strand_togetherness==r)); % I want this to sum faults that are going to be considered together
%fault_width_for_smoothing(r,1) = sum(fault_width(fault_strand_togetherness==r));
fault_width_for_smoothing(r,1) = fault_width(r,1); % assuming that all faults are the same depth
first_patch_in_strand_for_smoothing_master(r,1) = find(fault_strand_identifyer_for_smoothing==r,1);
last_patch_in_strand_for_smoothing_master(r,1) = find(fault_strand_identifyer_for_smoothing==r,1, 'last');
n_down_dip_patches_for_smoothing(r,1) = total_n_down_dip_patches; % because we're treating it as one fault. assuming all faults have some number of down dip patches
n_along_strike_patches_for_smoothing(r,1) = sum(n_along_strike_patches(fault_strand_togetherness==r)); % because we're treating it as one fault
n_slip_patches_on_each_fault_strand_for_smoothing(r,1) = sum(fault_strand_identifyer_for_smoothing==r);
end
elseif strcmp(invert.smooth_across_fault_strands, 'no') == 1
n_fault_strands_for_smoothing = n_fault_strands;
n_down_dip_patches_for_smoothing = n_down_dip_patches; % because we're treating it as separate faults as described in fault descriptor file
n_along_strike_patches_for_smoothing = n_along_strike_patches; % because we're treating it as separate faults as described in fault descriptor file
first_patch_in_strand_for_smoothing_master = first_patch_in_strand;
last_patch_in_strand_for_smoothing_master = last_patch_in_strand;
fault_length_for_smoothing = fault_length; % this is one value for each separate strand, regardless of whether attached or not
fault_width_for_smoothing = fault_width; % this is one value for each separate strand, regardless of whether attached or not
n_slip_patches_on_each_fault_strand_for_smoothing = n_slip_patches_on_each_fault_strand;
end
first_patch_in_strand_for_smoothing = first_patch_in_strand_for_smoothing_master;
last_patch_in_strand_for_smoothing = last_patch_in_strand_for_smoothing_master;
n_slip_patches_ON_on_each_fault_strand_for_smoothing = n_slip_patches_on_each_fault_strand_for_smoothing; % because they're all on to start off with
previous_n_slip_patches_ON_on_each_fault_strand_for_smoothing = n_slip_patches_ON_on_each_fault_strand_for_smoothing;
% *********** % ************ % ********** % ********** % ********** % *****
% If in testing mode...
% *********** % ************ % ********** % ********** % ********** % *****
% Sort out synthetic data...................................
if strcmp(testing.testing_mode, 'yes') == 1
% synthetic data.......................................................
load(testing.making_model, 'u'); % these are the true observations, from the forward model of the true slip
d_E = u(1,:);
d_N = u(2,:);
d_up = u(3,:);
d = [d_E'; d_N'; d_up']; % THIS IS d in E, N, up , as opposed to LOS
n_data = length(d);
d_identifyer = zeros(length(d), 1);
% add noise, if you like...............................................
if testing.add_noise ~= 0 % add_noise tells you the std of noise to add. add_noise == 1 means add 1 sigma noise, add_noise == 2 means add 2 sigma noise
% calculate standard deviation
sigmas = std(u');
sigma_E = sigmas(1);
sigma_N = sigmas(2);
sigma_up = sigmas(3);
% create an array of random numbers, in the sigma range. multiply by add_noise, so that you add the right number of sigmas
noise_E = (-sigma_E + (sigma_E - -sigma_E) * rand( length(u),1)) * testing.add_noise;% create an array of random numbers, in the sigma range
noise_N = (-sigma_N + (sigma_N - -sigma_N) * rand( length(u),1)) * testing.add_noise;% create an array of random numbers, in the sigma range
noise_up = (-sigma_up + (sigma_up - -sigma_up) * rand( length(u),1)) * testing.add_noise;% create an array of random numbers, in the sigma range
% update data vector
d_E = u(1,:) + noise_E';
d_N = u(2,:) + noise_N';
d_up = u(3,:) + noise_up';
d = [d_E'; d_N'; d_up'];
end
d_E = []; % get rid of it so you don't get confused / mess it up with removing point a mindist from the fault
d_N = [];
d_up = [];
% find out true moment.................................................
true1 = load(testing.making_model, 'spatial_model2');
true2 = load(testing.making_model, 'spatial_model3');
true3 = load(testing.making_model, 'synthetic_slip');
true3.synthetic_slip = reshape(true3.synthetic_slip, length(true3.synthetic_slip), 1);
true4 = load(testing.making_model, 'total_n_slip_patches');
M0_true = sum((elastic_params.mu_okada) * (reshape(true1.spatial_model2, true4.total_n_slip_patches,1) .* reshape(true2.spatial_model3, true4.total_n_slip_patches,1)) .* (true3.synthetic_slip));
data.seismic_moment = M0_true;
end
% *********** % ************ % ********** % ********** % ********** % *****
% Create the G kernel............................................
% *********** % ************ % ********** % ********** % ********** % *****
disp('... and calculating G, using Okada 1985...')
if strcmp(testing.testing_mode, 'yes') == 1 % OR IF ONLY USING GPS DATA, basically this is just not projected into LOS
for i = 1:total_n_slip_patches
[u_for_unit_slip, ~]=disloc3d3(disloc_model(:,i), locs, elastic_params.lambda,elastic_params.mu_okada) ;
G_E(:,i) = u_for_unit_slip(1,:);
G_N(:,i) = u_for_unit_slip(2,:);
G_up(:,i) = u_for_unit_slip(3,:);
G = [G_E; G_N; G_up];
end
end
if strcmp(invert.solve_for_dip, 'yes') == 1
dip_LUT = priors.min_dip : priors.max_dip; % Make a dip look up table
elseif strcmp(invert.solve_for_dip, 'no') == 1
dip_LUT = disloc_model(4);
end
if strcmp(invert.variable_rake, 'yes') == 1
G_ss = zeros( n_data, total_n_slip_patches, length(dip_LUT));
G_ds = zeros( n_data, total_n_slip_patches, length(dip_LUT));
for i = 1: length(dip_LUT)
% Calculate G for purely strike-slip
disloc_model_ss = disloc_model;
disloc_model_ss(5,:) = 0; % set rake of every patch to 0, which is purely LEFT LATERAL STRIKE SLIP
disloc_model_ss(4,:) = dip_LUT(i); % look up for each value of dip in LUT
%[ G_ss(:,:,i), ~, ~, ~ ] = calculate_G( total_n_slip_patches, testing, data, disloc_model_ss, locs_struct, elastic_params, los_vector_struct );
for j = 1: total_n_slip_patches
[u_for_unit_slip, ~]=disloc3d3(disloc_model_ss(:,j), locs, elastic_params.lambda, elastic_params.mu_okada) ;
if strcmp(testing.testing_mode, 'yes') ==1
G_E_ss(:,j) = u_for_unit_slip(1,:);
G_N_ss(:,j) = u_for_unit_slip(2,:);
G_up_ss(:,j) = u_for_unit_slip(3,:);
G_ss = [G_E_ss; G_N_ss; G_up_ss];
else
G_ss(:,j) = sum(u_for_unit_slip.*los_vector);
end
end
% Calculate G for purely dip-slip
disloc_model_ds = disloc_model;
disloc_model_ds(5,:) = 90; % set rake of every patch to 90, which is purely THRUST
disloc_model_ds(4,:) = dip_LUT(i);
%[ G_ds(:,:,i), ~, ~, ~] = calculate_G( total_n_slip_patches, testing, data, disloc_model_ds, locs_struct, elastic_params, los_vector_struct );
for j = 1: total_n_slip_patches
[u_for_unit_slip, ~]=disloc3d3(disloc_model_ds(:,j), locs, elastic_params.lambda, elastic_params.mu_okada) ;
if strcmp(testing.testing_mode, 'yes') ==1
G_E_ds(:,j) = u_for_unit_slip(1,:);
G_N_ds(:,j) = u_for_unit_slip(2,:);
G_up_ds(:,j) = u_for_unit_slip(3,:);
G_ds = [G_E_ds; G_N_ds; G_up_ds];
else
G_ds(:,j) = sum(u_for_unit_slip.*los_vector);
end
end
end
else
%[ G, G_E, G_N, G_up ] = calculate_G( total_n_slip_patches, testing, data, disloc_model, locs_struct, elastic_params, los_vector_struct );
for i = 1: total_n_slip_patches
[u_for_unit_slip, ~]=disloc3d3(disloc_model(:,i), locs, elastic_params.lambda, elastic_params.mu_okada) ;
G(:,i) = sum(u_for_unit_slip.*los_vector);
end
G_ss = [];
G_ds = [];
end
% *********** % ************ % ********** % ********** % ********** % *****
% Variance-covariance matrix
% *********** % ************ % ********** % ********** % ********** % *****
% Set up data hyperparameter
beta_step_size_initial = [];
beta_ratio = 1;
if strcmp(invert.solve_for_beta, 'yes') == 1
beta_initial = invert.beta_initial;
beta_step_size_initial = invert.beta_step_size;
beta_keep = zeros(1,ceil(invert.iterations));
max_beta = priors.max_beta;
min_beta = priors.min_beta;
else
max_beta = [];
min_beta = [];
end
if strcmp(testing.testing_mode, 'no') == 1
%if you have InSAR data ..................................................
if strcmp(data.InSAR_datafile, 'none') ~= 1
disp(' ');
disp(['Calculating var-covar matrix for Insar file ', data.InSAR_datafile]);
% Assuming previously calculated sill and nugget and range, so load it.
sigma_d_InSAR = [];
for p = 1:size(data.varcovar_details,2)
[sill, nugget, variogram_range]=textread(cell2mat(data.varcovar_details(p)),'%f %f %f'); % range in meters
[npoints]=textread(cell2mat(data.quadtree_n_points(p)),'%f'); % range in meters
% Calculate distance between all points
x = locs_InSAR(1,InSAR_identifyer==p); % meters
y = locs_InSAR(2,InSAR_identifyer==p);
x_dists = pdist2(x',x');
y_dists = pdist2(y',y');
dist = sqrt( x_dists.^2 + y_dists.^2); % this is in meters
% Calculate var-covar matrix
sigma_d_tmp = (sill-nugget) * exp((-3*dist)/variogram_range)+diag(nugget./ npoints); %if using quadtree.
sigma_d_InSAR = blkdiag(sigma_d_InSAR, sigma_d_tmp);
end
% If you you want to weight it ....................................................
sigma_d_InSAR = sigma_d_InSAR/data.weight_InSAR;
else
sigma_d_InSAR = [];
end
% if you have GPS data and you want to weight it ....................................................
if strcmp(data.GPS_datafile, 'none') ~= 1
sigma_d_GPS = diag(sigma_gps.^2) / data.weight_GPS; % variance is standard deviation squared. so if InSAR weight = 1 and GPS weight = 2, then we half the errors on GPS, coz we trust it quite as much
else
sigma_d_GPS = [];
end
% if you have atolls data and you want to weight it ................................................
if strcmp(data.atolls_datafile, 'none') ~= 1
sigma_d_atolls = diag(sigma_atolls.^2) / data.weight_atolls; % so if InSAR weight = 1 and GPS weight = 2, then we half the errors on GPS, coz we trust it quite as much
else
sigma_d_atolls = [];
end
sigma_d = blkdiag(sigma_d_InSAR, sigma_d_GPS, sigma_d_atolls);
elseif strcmp(testing.testing_mode, 'yes') == 1
datastd = testing.std * ones(length(d),1);
sigma_d = diag(datastd).^2; % assume no covar, only diagonals
end
inv_sigma_d = inv(sigma_d);
%% Invert for slip - it's why we're all here, this is the good bit ********
disp(' ');
disp('Solving for slip...');
disp('Oh boy, this is going to be fun.');
disp(' ');
%% OPTION 1 - least-squares inversion *************************************
% AT THE MOMENT this doesn't solve for variable rake... G is calculated
% just using the rake given in fault.fault_descriptor_file
if strcmp(invert.inversion_type, 'least_squares') ==1 && strcmp(invert.smoothing, 'laplacian') == 1
disp('Finding Laplacian matrix, with smoothing factor = ')
[smooth_los_kernel,smooth_obsdispl_los] = add_smoothing(G, d, spatial_model1, spatial_model2, spatial_model3, spatial_model4, (invert.scalar_smoothing_factor^2), n_datasets) ; % borrowed from gjf. This function takes an already existing kernel and observation set, and adds a Laplacian-minimising smoothing operation to it.
G=smooth_los_kernel ;
d=smooth_obsdispl_los ;
end
if strcmp(invert.inversion_type, 'least_squares') == 1 && strcmp(invert.smoothing, 'tikhonov') + strcmp(invert.smoothing, 'cv') ~= 1
disp('Solving for slip using fast non-negative least-squares...');
ATA= G'* G;
ATD= G'*d ;
[least_sq_solution, ~]=fnnls(ATA,ATD) ;
least_sq_residuals = calc_WRMS_offset( least_sq_solution, G, d, inv_sigma_d, offset_curr );
end
% For Tikhonov solution we want to minimise
% (d-Gm)^2 + mu*(m)^2
% Which is the damped least squares solution with the special case of
% minimising the weighted misfit of hte data misfit and the model norm. This
% is the Tikhonov solution. To find the value of mu we solve the solution
% repeatedly till we find the value such that
% (d-Gm)^2 = N * sigma^2
% (d-Gm) = sqrt(N) * sigma
% Where N is the number of degrees of freedom and sigma^2 is the variance
% (Since on average each error should equal sigma^2. So for N measurements
% the error should be N*sigma^2.)
% The solution for the damped least-squares solution is:
% m_tik = inv(G'*G + mu*I) * (G'*d)
% Or if you want to include data error:
% m_tik = inv(G' * inv_sigma_d * G + mu*I) * (G' * inv_sigma_d * d)
% if strcmp(invert.inversion_type, 'least_squares') == 1 && strcmp(invert.smoothing, 'tikhonov') == 1
% disp('Solving for slip using a linear least-squares tikhonov inversion...');
% %G_tik = [G, ones(length(d),1) ]; % coz I want to solve for offset too
% G_tik = G; % coz I want to solve for offset too
% mew_tests = [0.1:0.1:10, 10:1:100];
% tikhonov_criteria = total_n_slip_patches * sigma_d(1);
% %tikhonov_criteria = sqrt(total_n_slip_patches) * sqrt(sigma_d(1));
%
% for i= 1:length(mew_tests)
% mew = mew_tests(i);
% m_tik(:,i) = (G_tik' * inv_sigma_d * G_tik + mew * eye(total_n_slip_patches,total_n_slip_patches)) \ G_tik' * inv_sigma_d * d; % since inv(G' * G + mu * I)*(G' * d) is slower and less accurate than (G' * G + mu * I)\(G' * d)
% resid(:,i) = sum((d - (G_tik*m_tik(:,i))).^2);
% %resid(:,i) = sum(d - (G_tik*m_tik(:,i)));
% tik_test(1,i) = resid(:,i)-tikhonov_criteria;
% end
%
% % The best result is the one where (d-Gm) is closest to sqrt(N)*sigma
% [~,I] = min(abs(resid-tikhonov_criteria));
% best_mew = mew_tests(I);
% best_m_tik = m_tik(:,I);
% best_slip_tik = best_m_tik(1:end);
%
% figure; plot( mew_tests, resid )
% ylabel('sum( d - Gm)')
% xlabel('mew')
%
% best_offset_tik = 0;
% end
%%%%%%%%%%%%%%% This is solving for rake as well %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% but no positivity constraint %%%%%%%%%%%%%%%%%%%%%%%%
% if strcmp(invert.inversion_type, 'least_squares') == 1 && strcmp(invert.smoothing, 'tikhonov') == 1
% disp('Solving for slip and rake using a linear least-squares tikhonov inversion...');
% G_tik = [G_ss, G_ds ];
% mew_tests = [0:10:1000];
% m_tik = zeros(total_n_slip_patches*2, length(mew_tests));
% resid = zeros(1, length(mew_tests));
% for i= 1:length(mew_tests)
% mew = mew_tests(i);
% m_tik(:,i) = (G_tik' * inv_sigma_d * G_tik + mew * eye(total_n_slip_patches*2,total_n_slip_patches*2)) \ (G_tik' * inv_sigma_d * d); % since inv(G' * G + mu * I)*(G' * d) is slower and less accurate than (G' * G + mu * I)\(G' * d)
% slip_tik = m_tik(:,i);
% %offset_tik = m_tik(end);
% least_sq_residuals = calc_WRMS_offset( slip_tik, G_tik, d, inv_sigma_d, zeros(length(d),1) );
% resid(:,i) = sum((d - G_tik*m_tik(:,i)).^2);
% %figure; faults= disloc_model; faults(6,:) = sqrt(m_tik(1:total_n_slip_patches,i).^2 + m_tik(total_n_slip_patches+1:end,i).^2); doplot3d(faults', 'jet');
% end
%
% %The best result is the one where (d-Gm) is closest to sqrt(N)*sigma
% tikhonov_criteria = (length(m_tik(:,1))) * sigma_d(1);
% [~,I] = min(abs(resid-tikhonov_criteria));
% best_mew = mew_tests(I);
% best_m_tik = m_tik(:,I);
% best_ss_slip_tik = best_m_tik(1:total_n_slip_patches);
% best_ds_slip_tik = best_m_tik(total_n_slip_patches+1:end);
% best_slip_tik = sqrt(best_ss_slip_tik.^2 + best_ds_slip_tik.^2);
% best_rake_tik = 180 + atand( best_ds_slip_tik./best_ss_slip_tik);
%
% theta = 180 - best_rake_tik;
% theta = degtorad(theta);
% for r = 1 : total_n_slip_patches % we wanna select the correct COLUMN, relating to the current slip patch, out of the G matrix for that value of rake.
% best_G_tik(:,r) = cos(theta(r)) * G_ss(:,r) + sin(theta(r)) * G_ds(:,r); % I THINK. check yo' trigonometry
% end
%
% best_offset_tik = 0;
% end
%%%%%%%%%%%%%%% This is tikhonov using fnnls %%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(invert.inversion_type, 'least_squares') == 1 && strcmp(invert.smoothing, 'tikhonov') == 1
disp('Solving for slip and rake using a linear least-squares tikhonov inversion...');
% We calculate Gss and Gds above, but I've copied this here, to
% solve for 45 degrees around a specified rake
centrerake = disloc_model(5,1);
minrake = centrerake-45;
maxrake = centrerake+45;
disloc_model_minrake = disloc_model;
disloc_model_minrake(5,:) = minrake; % set rake of every patch to 0, which is purely RIGHT LATERAL STRIKE SLIP
disloc_model_minrake(4,:) = dip_LUT(i); % look up for each value of dip in LUT
[ G_minrake(:,:,i), ~, ~, ~ ] = calculate_G( total_n_slip_patches, testing, data, disloc_model_minrake, locs_struct, elastic_params, los_vector_struct );
% [ G_ss(first_patch_in_strand(i):last_patch_in_strand(i),first_patch_in_strand(i):last_patch_in_strand(i)), ~, ~, ~ ] = calculate_G( total_n_slip_patches, testing, data, disloc_model_ss, locs_struct, elastic_params, los_vector_struct );
% Calculate G for purely dip-slip
disloc_model_maxrake = disloc_model;
disloc_model_maxrake(5,:) = maxrake; % set rake of every patch to 90, which is purely THRUST
disloc_model_maxrake(4,:) = dip_LUT(i);
[ G_maxrake(:,:,i), ~, ~, ~] = calculate_G( total_n_slip_patches, testing, data, disloc_model_maxrake, locs_struct, elastic_params, los_vector_struct );
mew_tests = 1:10:100;
m_tik = zeros(total_n_slip_patches*2, length(mew_tests));
resid = zeros(1, length(mew_tests));
d_tik = [d; zeros((total_n_slip_patches*2),1)];
weights = chol(inv(sigma_d));
for i= 1:length(mew_tests)
% Make new G using smoothing factor
mew = mew_tests(i);
G = [G_minrake, G_maxrake ; mew * eye(total_n_slip_patches*2) ];
% Add weights and apply
allW_tik = blkdiag(weights,eye(total_n_slip_patches*2)); % add in ones along diagonal for minimum norm
G_tik = allW_tik*G;
obs_tik = allW_tik*d_tik;
% Solve for slip
ATA=G_tik'*G_tik ;
ATD=G_tik'*obs_tik ;
[m_tik(:,i), ~]=fnnls(ATA,ATD) ;
least_sq_residuals = calc_WRMS_offset( m_tik(:,i), [G_minrake, G_maxrake], d, inv_sigma_d, zeros(length(d),1) );
%resid(:,i) = sum((d - G_tik*m_tik(:,i)).^2);
%figure; faults= disloc_model; faults(6,:) = sqrt(m_tik(1:total_n_slip_patches,i).^2 + m_tik(total_n_slip_patches+1:end,i).^2); doplot3d(faults', 'jet');
end
%The best result is the one where (d-Gm) is closest to sqrt(N)*sigma
tikhonov_criteria = (length(m_tik(:,1))) * sigma_d(1);
[~,I] = min(abs(least_sq_residuals-tikhonov_criteria));
best_mew = mew_tests(I);
best_m_tik = m_tik(:,I);
best_minrake_slip_tik = best_m_tik(1:total_n_slip_patches);
best_maxrake_slip_tik = best_m_tik(total_n_slip_patches+1:end);
best_slip_tik = sqrt(best_minrake_slip_tik.^2 + best_maxrake_slip_tik.^2);
best_rake_tik = minrake + atand( best_maxrake_slip_tik./best_minrake_slip_tik);
best_rake_tik(best_minrake_slip_tik==0 & best_maxrake_slip_tik==0) = centrerake;
% theta = 180 - best_rake_tik;
% theta = degtorad(theta);
% for r = 1 : total_n_slip_patches % we wanna select the correct COLUMN, relating to the current slip patch, out of the G matrix for that value of rake.
% best_G_tik(:,r) = cos(theta(r)) * G_ss(:,r) + sin(theta(r)) * G_ds(:,r); % I THINK. check yo' trigonometry
% end
best_G_tik = G_ss*diag(cosd(best_rake_tik)) + G_ds*diag(sind(best_rake_tik));
best_offset_tik = 0;
end
%%%%%%%%%%%%%%%% SVD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Explanation of SVD goes here
if strcmp(invert.inversion_type, 'SVD') == 1 == 1
disp('Solving for slip and rake using SVD...');
G_svd = [G_ss, G_ds ];
[U, S, V] = svd(G_svd);
p = sum(diag(S)<1e-6);
Sp = S(1:p, 1:p);
Vp = V(:, 1:p);
Up = U(:, 1:p);
Gg = Vp * inv(Sp) * Up.';
m_svd = Gg * d;
least_sq_residuals = calc_WRMS_offset( m_svd, G_svd, d, inv_sigma_d, zeros(length(d),1) );
ss_slip_svd = m_svd(1:total_n_slip_patches);
ds_slip_svd = m_svd(total_n_slip_patches+1:end);
slip_svd = sqrt(ss_slip_svd.^2 + ds_slip_svd.^2);
rake_svd = 180 + atand(ds_slip_svd./ss_slip_svd);
% theta = 180 - rake_svd;
% theta = degtorad(theta);
% for r = 1 : total_n_slip_patches % we wanna select the correct COLUMN, relating to the current slip patch, out of the G matrix for that value of rake.
% best_G_svd(:,r) = cos(theta(r)) * G_ss(:,r) + sin(theta(r)) * G_ds(:,r); % I THINK. check yo' trigonometry
% end
best_G_svd = G_ss*diag(cosd(rake_svd)) + G_ds*diag(sind(rake_svd));
best_offset_tik = 0;
end
%%%%%%%%%%%%%%%%% Cross validation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For cross-validation we pick which is the best smoothing factor by
% measuring the fit to missing data. A data point is dropped, the minimum
% norm solution is calculated for a given range of smoothing factors, and
% then the best smoothing factor is the one that predicts the missing data
% best.
% Tom has altered this (and sped it up) by instead of dropping one data
% point each time, he instead randomly reorders the data vector
% and cuts into k subsets. He uses k-1 subsets to estimate the inversion
% parameters and predict the data values at the remaining subset. This is
% done k times for each smoothing factor and the total misfit calculated by