doc updates for SeisInvExpt

docs/Manifest.toml

@@ -276,6 +276,12 @@ version = "7.0.1"
 deps = ["Libdl"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
+[[LinearMaps]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "ecd59376591724b29db3728f54ee885ecbeb2395"
+uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
+version = "2.5.2"
+
 [[Loess]]
 deps = ["Distances", "Random", "Statistics", "Test"]
 git-tree-sha1 = "0ee46caf683a422b595be4dfaed6cda28f541e25"
@@ -367,9 +373,9 @@ uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[PlotThemes]]
 deps = ["PlotUtils", "Requires", "Statistics"]
-git-tree-sha1 = "5bf458e66ca6d617e299a55dfb5328bb26dbf731"
+git-tree-sha1 = "df772cc7c78862da96af1ee85cd0111c6640e44e"
 uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
-version = "1.0.0"
+version = "1.0.1"
 
 [[PlotUtils]]
 deps = ["Colors", "Dates", "Printf", "Random", "Reexport"]
@@ -493,9 +499,9 @@ version = "0.32.0"
 
 [[StatsFuns]]
 deps = ["Rmath", "SpecialFunctions"]
-git-tree-sha1 = "938f2520e6c093aaf5a7e7c38489e03c27a8c136"
+git-tree-sha1 = "79982835d2ff3970685cb704500909c94189bde9"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "0.9.2"
+version = "0.9.3"
 
 [[SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "SparseArrays"]

docs/Project.toml

@@ -1,5 +1,9 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/make.jl

@@ -1,8 +1,9 @@
 using Pkg; Pkg.add("Gadfly"); Pkg.add("Plots");
 
-using Documenter, GeoPhyInv
+using Documenter, GeoPhyInv, SparseArrays, LinearAlgebra, Random, LinearMaps
 using Gadfly, Plots
 gr()
+theme(:juno)
 
 #makedocs(modules=[GeoPhyInv],
 #	         doctest=true)
@@ -23,12 +24,16 @@ makedocs(
 				  "Tutorial I" => "generated/fdtd/create_snaps.md",
 				  "Tutorial II" => "generated/fdtd/reuse_expt.md",
 				  ],
+	    "SeisInvExpt" => Any[
+				  "Description" => "generated/fwi/doc.md",
+				  "Tutorial" => "generated/fwi/pizza.md",
+				  "Adjoint Tests" => "generated/fwi/born_tutorial.md",
+				 ],
 	    "PoissonExpt" => Any[
-				  "Introduction" => "generated/Poisson/doc.md",
-				 "Record data" => "generated/Poisson/forw.md",
-				 "Born map" => "generated/Poisson/test_born.md",
+				  "Description" => "generated/Poisson/doc.md",
+				 "Tutorial" => "generated/Poisson/forw.md",
+				 "Adjoint Tests" => "generated/Poisson/test_born.md",
 				 ],
-#	    "SeisInvExpt" => Any[]
 	    ]
 #	    "Seismic Born Modeling" => "FWI/born_map.md",
 #	    "Seismic Full Waveform Inversion" => "FWI/gradient_accuracy.md",

docs/src/generated/fdtd/reuse_expt.md

@@ -82,22 +82,22 @@ nothing #hide
 ```@example reuse_expt
 @time update!(pa);
 d1=copy(pa[:data][:P])
-p1=heatmap(pa[:data]);
+p1=heatmap(pa[:data], grid=true);
 nothing #hide
 ```
 
-### Update medium in `pa` without memory allocation
+### Change medium in `pa` without memory allocation
 
 ```@example reuse_expt
 update!(pa, medium_new)
 ```
 
-### Modeling #2 and plot
+### Modeling #2
 
 ```@example reuse_expt
 @time update!(pa);
 d2=copy(pa[:data][:P])
-p2=heatmap(pa[:data]);
+p2=heatmap(pa[:data], grid=true);
 nothing #hide
 ```
 
@@ -107,7 +107,7 @@ Test
 @test d1≠d2
 ```
 
-### Plotting
+### Plotting #2
 
 ```@example reuse_expt
 plot(p1,p2, size=(500, 300))

docs/src/generated/fwi/born_tutorial.md

@@ -0,0 +1,62 @@
+```@meta
+EditURL = "<unknown>/fwi/born_tutorial.jl"
+```
+
+```@example born_tutorial
+using Test
+using GeoPhyInv
+using LinearMaps
+using Misfits
+using LinearAlgebra
+```
+
+### SeisInvExpt
+
+```@example born_tutorial
+pa=SeisInvExpt(FdtdBorn(), LS(), :pizza, rfields=[:P, :Vz]); # use linearized modeling attribute
+F=LinearMap(pa); # generate LinearMap
+nothing #hide
+```
+
+### Linearity test
+
+```@example born_tutorial
+x1=randn(size(F,2)); # random input 1
+x2=randn(size(F,2)); # random input 2
+x12=x1.+x2; # sum
+nothing #hide
+```
+
+```@example born_tutorial
+d12=F*x12;
+d1=F*x1; #
+d2=F*x2;
+d12new=d1.+d2;
+nothing #hide
+```
+
+```@example born_tutorial
+@test Misfits.error_squared_euclidean!(nothing, d12, d12new, nothing, norm_flag=true)<1e-25
+```
+
+### Adjoint test
+
+```@example born_tutorial
+x=randn(size(F,2)); # random x
+y=randn(size(F,1)); # random y
+a=LinearAlgebra.dot(y,F*x);
+b=LinearAlgebra.dot(x,adjoint(F)*y);
+c=LinearAlgebra.dot(x, transpose(F)*F*x);
+nothing #hide
+```
+
+```@example born_tutorial
+println("adjoint test: ", a, "\t", b)
+@test isapprox(a,b,rtol=1e-5)
+```
+
+```@example born_tutorial
+println("x*transpose(F)*F*x be positive: ", c)
+@test c>0.0
+```
+

docs/src/generated/fwi/doc.md

@@ -0,0 +1,23 @@
+```@meta
+EditURL = "<unknown>/fwi/doc.jl"
+```
+
+```@example doc
+using BenchmarkTools
+using GeoPhyInv
+using Test
+```
+
+# Intro
+
+```@docs
+SeisInvExpt
+Base.getindex(::GeoPhyInv.PFdtd, ::Symbol)
+```
+
+# Methods
+```@docs
+update!(::GeoPhyInv.PFWI)
+LinearMaps.LinearMap(::GeoPhyInv.PFWI)
+```
+

docs/src/generated/fwi/pizza.md

@@ -0,0 +1,54 @@
+```@meta
+EditURL = "<unknown>/fwi/pizza.jl"
+```
+
+```@example pizza
+using GeoPhyInv
+using Plots; gr();
+nothing #hide
+```
+
+### Medium
+
+```@example pizza
+medium_true=Medium(:pizza)
+```
+
+### SeisInvExpt
+
+```@example pizza
+pa=SeisInvExpt(Fdtd(), LS(), :pizza);
+nothing #hide
+```
+
+### Initial plotting
+
+```@example pizza
+p1=heatmap(pa.mediumm, :vp, title="Initial vp")
+scatter!(pa.ageom, SSrcs())
+scatter!(pa.ageom, Recs())
+p2=heatmap(medium_true, :vp, title="True vp")
+scatter!(pa.ageom, SSrcs())
+scatter!(pa.ageom, Recs())
+plot(p2, p1, size=(700,300))
+```
+
+### Final plotting
+
+```@example pizza
+update!(pa,solver=:ipopt);
+nothing #hide
+```
+
+### Plotting #1
+
+```@example pizza
+p1=heatmap(pa.mediumm, :vp, title="Inverted vp")
+scatter!(pa.ageom, SSrcs())
+scatter!(pa.ageom, Recs())
+p2=heatmap(medium_true, :vp, title="True vp")
+scatter!(pa.ageom, SSrcs())
+scatter!(pa.ageom, Recs())
+plot(p2, p1, size=(700,300))
+```
+

docs/src/generated/media/doc.md

@@ -44,6 +44,13 @@ Inspect the bounds.
 mod.bounds
 ```
 
+Plotting is easy
+
+```@example doc
+p1=heatmap(mod, :vp)
+plot(p1,size=(800,400))
+```
+
 To construct an instance of `Medium`, we need a 2-D grid.
 
 ```@example doc
@@ -95,18 +102,6 @@ update!(mod, [:vp,:rho], randn_perc=1.);
 nothing #hide
 ```
 
-Some plotting #1
-
-```@example doc
-p1=heatmap(mod, :vp); plot(p1)
-```
-
-Some plotting #2
-
-```@example doc
-p2=heatmap(mod, :vs); plot(p2)
-```
-
 # Methods
 ```@docs
 Base.getindex(::Medium, ::Symbol)

docs/src/generated/srcwav/doc.md

@@ -74,11 +74,11 @@ nothing #hide
 Most of the methods listed below are also applicable to individual elements of `srcwav`.
 
 ```@docs
-update!(GeoPhyInv.VNamedD, ::Vector{Symbol}, ::AbstractVector{Float64},)
+update!(GeoPhyInv.VNamedD, ::Vector{Symbol}, ::AbstractArray)
 Base.reverse!(::GeoPhyInv.VNamedD)
 Base.iszero(::GeoPhyInv.VNamedD)
 Base.isequal(::GeoPhyInv.VNamedD, ::GeoPhyInv.VNamedD)
-issimilar(::GeoPhyInv.VNamedD, ::GeoPhyInv.VNamedD)
+GeoPhyInv.issimilar(::GeoPhyInv.VNamedD, ::GeoPhyInv.VNamedD)
 Base.vec(::GeoPhyInv.VNamedD)
 Random.randn!(::GeoPhyInv.VNamedD)
 Base.fill!(::GeoPhyInv.VNamedD, ::Float64)

src/database/core.jl

@@ -103,7 +103,7 @@ update!(srcwav[1], [:P, :Vx], w)
 ```
 Populate first supersource of `srcwav` by a wavelet `w` for `:P` and `:Vx` fields.
 """
-function update!(d::NamedD, fields::Vector{Symbol}, w::AbstractVector{Float64},)
+function update!(d::NamedD, fields::Vector{Symbol}, w::AbstractArray)
 	@assert length(w)==length(d.grid)
 	@inbounds for f in fields
 		@assert f ∈ names(d.d)[1]
@@ -116,13 +116,15 @@ function update!(d::NamedD, fields::Vector{Symbol}, w::AbstractVector{Float64},)
 	return d
 end
 
+
+
 """
 ```julia
 update!(srcwav, [:P, :Vx], w)
 ```
 Populate `srcwav` by a wavelet vector `w` for `:P` and `:Vx` fields.
 """
-function update!(dat::VNamedD, fields::Vector{Symbol}, w::AbstractVector{Float64},)
+function update!(dat::VNamedD, fields::Vector{Symbol}, w::AbstractArray)
 	for d in dat
 		update!(d,fields,w)
 	end
@@ -155,6 +157,7 @@ Return if two `NamedD`'s have same dimensions and bounds.
 function issimilar(dat1::NamedD, dat2::NamedD)
 	return isequal(dat1.grid, dat2.grid) && isequal(names(dat1.d), names(dat2.d)) && (size(dat1.d)==size(dat2.d)) 
 end
+
 """
 ```julia
 issimilar(srcwav1, srcwav2)
@@ -354,7 +357,6 @@ fill!(srcwav, b)
 ```
 In-place method to fill `srcwav` with scalar `b`.
 """
-
 function Base.fill!(data::VNamedD, k::Float64)
 	for d in data
 		fill!(d,k)
@@ -376,7 +378,6 @@ reverse!(srcwav)
 ```
 Perform in-place time-reversal operation for each wavelet in `srcwav`.
 """
-
 function Base.reverse!(data::VNamedD)
 	for d in data
 		reverse!(d)

src/database/misfits.jl

@@ -1,8 +1,8 @@
 mutable struct VNamedD_misfit{T}
 	x::Vector{NamedD{T}} # modelled data
-	y::Vector{NamedD{T}} # modelled data
-	w::Vector{NamedD{T}} # modelled data
-	dJx::Vector{NamedD{T}} # modelled data
+	y::Vector{NamedD{T}} # observed data
+	w::Vector{NamedD{T}} # weights
+	dJx::Vector{NamedD{T}} # gradient 
 	ynorm::Float64 # normalize functional with energy of y
 end
 

src/fdtd/fdtd.jl

@@ -191,7 +191,7 @@ end
 
 """
 ```julia
-pa = SeisForwExpt(attrib_mod, ageom=ageom, srcwav=srcwav, medium=mod, tgrid=tgrid);
+pa = SeisForwExpt(attrib_mod; ageom, srcwav, medium, tgrid);
 ```
 Method to create an instance of `SeisForwExpt`. 
 The output of this method can be used as an input to the in-place method `update!`, to actually perform a
@@ -206,6 +206,9 @@ finite-difference modeling.
 * `tgrid` : modeling time grid, maximum time in `tgrid`should be greater than or equal to maximum source time, same sampling interval as in `srcwav`
 * `ageom` :  acquisition 
 * `srcwav` : source wavelets  
+
+# Optional Keyword Arguments 
+
 * `sflags=2` : source related flags 
   * `=0` inactive sources
   * `=1` sources with injection rate