Merge pull request #87 from slimgroup/scalar-reduce

Misc bug fixes

.github/workflows/cancel.yml

@@ -12,5 +12,5 @@ jobs:
         with:
           # Ids to cancel op/judi
           # https://api.github.com/repos/slimgroup/JUDI.jl/actions/workflows
-          workflow_id: 1223553, 1223567
+          workflow_id: 1223553, 1223567, 18418460
           access_token: ${{ github.token }}

.github/workflows/ci-examples.yml

@@ -0,0 +1,62 @@
+name: CI-examples
+
+on: 
+  # Trigger the workflow on push to master or pull request
+  # to be merged in master
+  push:
+    branches:
+      - master
+  pull_request:
+    branches:
+      - master
+
+jobs:
+  test:
+    name: JUDI example on Julia ${{ matrix.version }}
+    runs-on: ${{ matrix.os }}
+    env:
+      DEVITO_ARCH: gcc-9
+      DEVITO_LANGUAGE: "openmp"
+      DEVITO_BACKEND: "core"
+      DEVITO_LOGGING: "ERROR"
+      OMP_NUM_THREADS: ${{ matrix.omp }}
+      NITER: 1
+
+    strategy:
+      fail-fast: false
+
+      matrix:
+        os: [ubuntu-latest]
+        omp: [2]
+        version: [1.6, 1.7]
+
+    steps:
+      - name: Checkout JUDI
+        uses: actions/checkout@v2
+
+      - name: Install GCC 9
+        if : runner.os == 'macOS'
+        run : brew install gcc@9
+
+      - name: Set julia python
+        run: |
+          PYTHON=$(which python3) julia -e 'using Pkg;Pkg.add("PyCall");Pkg.build("PyCall")'
+
+      - name: Setup julia
+        uses: julia-actions/setup-julia@v1
+        with:
+          version: "${{ matrix.version }}"
+          arch: x64
+
+      - name: Build JUDI
+        uses: julia-actions/julia-buildpkg@latest
+
+      - name: Install packages
+        run: |
+          julia -e 'using Pkg;Pkg.add(["NLopt", "JOLI", "PyPlot", "IterativeSolvers", "SlimOptim", "HDF5", "SegyIO", "SetIntersectionProjection"])'
+          julia -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd()))'
+
+      - name: Run base examples
+        working-directory: examples/scripts
+        run: |
+          julia -p 2 -e 'for f in filter!(e->e≠"runall.jl", readdir(".")); include(f); end'

.github/workflows/ci-judi.yml

@@ -12,7 +12,7 @@ on:
 
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }}
+    name: JUDI base on Julia ${{ matrix.version }} - ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     env:
       DEVITO_ARCH: gcc-7
@@ -25,14 +25,8 @@ jobs:
       fail-fast: false
 
       matrix:
-        version: ['1.2', '1.3', '1.4', '1.5', '1.6', '1.7']
-        os: [ubuntu-latest]
-        include:
-          - version: 1.3
-            os: macos-latest
-
-          - version: 1.4
-            os: macos-latest
+        version: ['1.6', '1.7']
+        os: [ubuntu-latest, macos-latest]
 
     steps:
       - name: Checkout JUDI

.github/workflows/ci-op.yml

@@ -28,29 +28,19 @@ jobs:
       matrix:
         os: [ubuntu-latest]
         op: ["ISO_OP", "ISO_OP_FS", "TTI_OP", "TTI_OP_FS", "BASICS"]
-        version: ['1.3', '1.4', '1.5']
+        version: ['1.6', '1.7']
         omp: [2]
 
         include:
           - os: macos-latest
-            version: '1.4'
+            version: '1.6'
             op: "BASICS"
             omp: 1
 
-          - os: ubuntu-latest
+          - os: macos-latest
             version: '1.7'
             op: "ISO_OP_FS"
-            omp: 2
-
-          - os: ubuntu-latest
-            version: '1.7'
-            op: "ISSUES"
-            omp: 2
-
-          - os: ubuntu-latest
-            version: '1.6'
-            op: "ISO_OP"
-            omp: 2
+            omp: 1
 
     steps:
       - name: Checkout JUDI

examples/scripts/fwi_example_2D.jl

@@ -11,9 +11,9 @@ n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0
 model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)
 
 # Bound constraints
-v0 = sqrt.(1f0 ./ model0.m)
-vmin = ones(Float32,model0.n) .* 1.3f0
-vmax = ones(Float32,model0.n) .* 6.5f0
+v0 = sqrt.(1f0 ./ m0)
+vmin = ones(Float32, model0.n) .* 1.3f0
+vmax = ones(Float32, model0.n) .* 6.5f0
 vmin[:,1:21] .= v0[:,1:21]   # keep water column fixed
 vmax[:,1:21] .= v0[:,1:21]
 
@@ -34,9 +34,9 @@ q = judiVector(src_geometry,wavelet)
 F0 = judiModeling(deepcopy(model0), src_geometry, d_obs.geometry)
 
 # Optimization parameters
-niterations = 10
+niterations = parse(Int, get(ENV, "NITER", "10"))
 batchsize = 16
-fhistory_SGD = zeros(Float32,niterations)
+fhistory_SGD = zeros(Float32, niterations)
 
 # Projection operator for bound constraints
 proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)]; dims=2),model0.n)
@@ -47,7 +47,7 @@ for j=1:niterations
 
     # get fwi objective function value and gradient
     i = randperm(d_obs.nsrc)[1:batchsize]
-    fval, gradient = fwi_objective(model0,q[i],d_obs[i])
+    fval, gradient = fwi_objective(model0, q[i], d_obs[i])
     p = -gradient/norm(gradient, Inf)
 
     println("FWI iteration no: ",j,"; function value: ",fval)
@@ -57,9 +57,10 @@ for j=1:niterations
     function ϕ(α)
         F0.model.m .= proj(model0.m .+ α * p)
         misfit = .5*norm(F0[i]*q[i] - d_obs[i])^2
+        @show α, misfit
         return misfit
     end
-    step, fval = ls(ϕ, 1f0, fval, dot(gradient, p))
+    step, fval = ls(ϕ, 1f-1, fval, dot(gradient, p))
 
     # Update model and bound projection
     model0.m .= proj(model0.m .+ step .* p)

examples/scripts/fwi_example_NLopt.jl

@@ -11,7 +11,7 @@ n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0
 model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)
 
 # Bound constraints
-v0 = sqrt.(1f0 ./ model0.m)
+v0 = sqrt.(1f0 ./m0)
 vmin = ones(Float32, model0.n) .* 1.3f0
 vmax = ones(Float32, model0.n) .* 6.5f0
 
@@ -59,5 +59,5 @@ end
 opt = Opt(:LD_LBFGS, prod(model0.n))
 lower_bounds!(opt, mmin); upper_bounds!(opt, mmax)
 min_objective!(opt, f!)
-maxeval!(opt, 10)
+maxeval!(opt, parse(Int, get(ENV, "NITER", "10")))
 (minf, minx, ret) = optimize(opt, vec(model0.m.data))

examples/scripts/fwi_example_constraints.jl

@@ -1,19 +1,20 @@
 # 2D FWI on Overthrust model using minConf library
-# Author: Philipp Witte, [email protected]
+# Authors: 
+# Philipp Witte, [email protected]
 # Date: December 2017
-#
+# Mathias Louboutin, [email protected]
+# Date: January 2022
 
-using Statistics, Random, Pkg
-using LinearAlgebra
-using JUDI, SlimOptim, HDF5, SegyIO, PyPlot, FFTW
+using Statistics, Random, LinearAlgebra
+using JUDI, SlimOptim, HDF5, SegyIO, PyPlot
 using SetIntersectionProjection
 
 # Load starting model
 n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0")
 model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)
 
 # Bound constraints
-v0 = sqrt.(1f0 ./ model0.m)
+v0 = sqrt.(1f0 ./ m0)
 
 # Load data
 block = segy_read("../../data/overthrust_shot_records.segy")
@@ -27,7 +28,7 @@ q = judiVector(src_geometry,wavelet)
 ############################### FWI ###########################################
 
 # Optimization parameters
-niterations = 10
+niterations = parse(Int, get(ENV, "NITER", "10"))
 batchsize = 10
 fhistory_SGD = zeros(Float32,niterations)
 
@@ -50,7 +51,6 @@ options.rho_ini=[1.0f0]
 
 set_zero_subnormals(true)
 BLAS.set_num_threads(2)
-FFTW.set_num_threads(2)
 options.parallel=false
 options.feasibility_only = false
 options.zero_ini_guess=true
@@ -89,7 +89,7 @@ options.rho_ini = ones(length(TD_OP))*10.0
 
 proj_intersection = x-> PARSDMM(x, AtA, TD_OP, set_Prop, P_sub, model0, options)  
 
-function prj(input)
+function proj(input)
     input = Float32.(input)
     (x,dummy1,dummy2,dymmy3) = proj_intersection(vec(input.data))
     return reshape(x, model0.n)
@@ -118,7 +118,7 @@ for j=1:niterations
     step, fval = ls(ϕ, 1f0, fval, dot(gradient, p))
 
     # Update model and bound projection
-    model0.m = proj(model0.m .+ step .* p)
+    model0.m .= proj(model0.m .+ step .* p)
 end
 
 figure(); imshow(sqrt.(1f0./adjoint(model0.m))); title("FWI with SPG")

examples/scripts/fwi_example_minConf.jl

@@ -11,7 +11,7 @@ n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0
 model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)
 
 # Bound constraints
-v0 = sqrt.(1f0 ./ model0.m)
+v0 = sqrt.(1f0 ./ m0)
 vmin = ones(Float32,model0.n) .* 1.3f0
 vmax = ones(Float32,model0.n) .* 6.5f0
 vmin[:,1:21] .= v0[:,1:21]   # keep water column fixed
@@ -34,7 +34,7 @@ q = judiVector(src_geometry,wavelet)
 
 
 # Optimization parameters
-fevals = 16
+fevals = parse(Int, get(ENV, "NITER", "10"))
 batchsize = 8
 
 # Objective function for minConf library
@@ -52,13 +52,13 @@ function objective_function(x)
 end
 
 # Bound projection
-ProjBound(x) = median([mmin x mmax]; dims=2)
+proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)]; dims=2),model0.n)
 
 # FWI with SPG
 options = spg_options(verbose=3, maxIter=fevals, memory=3)
-sol = spg(objective_function, vec(model0.m), ProjBound, options)
+sol = spg(objective_function, vec(model0.m), proj, options)
 
 # Plot result
-imshow(reshape(sqrt.(1f0 ./ sol.sol), model0.n)', extent=[0, 10, 3, 0])
+imshow(reshape(sqrt.(1f0 ./ sol.x), model0.n)', extent=[0, 10, 3, 0])
 xlabel("Lateral position [km]")
 ylabel("Depth [km]")

examples/scripts/fwi_gauss_newton.jl

@@ -11,7 +11,7 @@ n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0
 model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)
 
 # Bound constraints
-v0 = sqrt.(1 ./ model0.m)
+v0 = sqrt.(1 ./ m0)
 vmin = ones(Float32,model0.n) .* 1.3f0
 vmax = ones(Float32,model0.n) .* 6.5f0
 vmin[:,1:21] .= v0[:,1:21]   # keep water column fixed
@@ -41,8 +41,8 @@ F = judiModeling(info,model0)
 J = judiJacobian(Pr*F*Ps',q)
 
 # Optimization parameters
-maxiter = 10
-maxiter_GN = 5
+maxiter = parse(Int, get(ENV, "NITER", "10"))
+maxiter_GN = parse(Int, get(ENV, "NITER", "5"))
 fhistory_GN = zeros(Float32,maxiter)
 proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)]; dims=2),model0.n)
 

examples/scripts/modeling_extended_source_3D.jl

@@ -30,7 +30,7 @@ model = Model(n, d, o, m)
 # Receiver geometry
 nxrec = 120
 nyrec = 100
-xrec = range(50f0, stop=1150f0, length=nxrec)
+xrec = range(50f0, stop=950f0, length=nxrec)
 yrec = range(100f0, stop=900f0, length=nyrec)
 zrec = 50f0
 
@@ -78,4 +78,4 @@ dw = adjoint(F)*d_sim
 # Jacobian
 J = judiJacobian(F, w)
 d_lin = J*dm
-g = adjoint(J)*d_lin
+g = adjoint(J)*d_lin

examples/scripts/modeling_medical_2D.jl

@@ -34,9 +34,9 @@ dtR = 0.2f0    # receiver sampling interval [ms]
 recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)
 
 ## Set up source geometry (cell array with source locations for each shot)
-xsrc = convertToCell([d[1]*61])
-ysrc = convertToCell([0f0])
-zsrc = convertToCell([d[1]])
+xsrc = d[1]*61
+ysrc = 0f0
+zsrc = d[1]
 
 # source sampling and number of time steps
 timeS = 250f0     # ms

src/TimeModeling/LinearOperators/judiLRWF.jl

@@ -81,7 +81,7 @@ adjoint(A::judiLRWF{DDT,RDT}) where {DDT,RDT} =
 # *(judiLRWF, judiWeights)
 function *(A::judiLRWF{ADDT,ARDT}, v::judiWeights{vDT}) where {ADDT,ARDT,vDT}
     A.n == size(v,1) || throw(judiLRWFexception("Shape mismatch: A:$(size(A)), v: $(size(v))"))
-    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiLRWF,judiVector):",A.name,typeof(A),vDT]," / "))
+    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiLRWF,judiWeights):",A.name,typeof(A),vDT]," / "))
     V = judiExtendedSource(A.info,A.wavelet,v.weights)
     jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiLRWF,judiWeights):",A.name,typeof(A),eltype(V)]," / "))
     return V
@@ -90,9 +90,9 @@ end
 # *(judiLRWF, vec)
 function *(A::judiLRWF{ADDT,ARDT}, v::AbstractVector{vDT}) where {ADDT,ARDT,vDT}
     A.n == size(v, 1) || throw(judiLRWFexception("Shape mismatch: A:$(size(A)), v: $(size(v))"))
-    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiLRWF,judiVector):",A.name,typeof(A),vDT]," / "))
+    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiLRWF,AbstractVector):",A.name,typeof(A),vDT]," / "))
     V = judiExtendedSource(A.info,A.wavelet, v)
-    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiLRWF,judiWeights):",A.name,typeof(A),eltype(V)]," / "))
+    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiLRWF,AbstractVector):",A.name,typeof(A),eltype(V)]," / "))
     return V
 end
 

src/TimeModeling/LinearOperators/judiPDE.jl

@@ -133,7 +133,7 @@ function *(A::judiPDE{ADDT,ARDT},v::judiRHS{vDT}) where {ADDT,ARDT,vDT}
     jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiRHS):",A.name,typeof(A),vDT]," / "))
 	args = (v.geometry, v.data, A.geometry, nothing, nothing)
     V = A.fop(args)
-    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
+    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiRHS):",A.name,typeof(A),eltype(V)]," / "))
     return V
 end
 
@@ -142,7 +142,7 @@ function *(A::judiPDEadjoint{ADDT,ARDT},v::judiRHS{vDT}) where {ADDT,ARDT,vDT}
     jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiRHS):",A.name,typeof(A),vDT]," / "))
 	args = (A.geometry,nothing,v.geometry,v.data,nothing)
     V = A.fop(args)
-    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
+    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiRHS):",A.name,typeof(A),eltype(V)]," / "))
     return V
 end
 
@@ -152,7 +152,7 @@ function *(A::judiPDE{ADDT,ARDT},v::judiWavefield{vDT}) where {ADDT,ARDT,vDT}
 	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiWavefield):",A.name,typeof(A),vDT]," / "))
 	args = (nothing,v.data,A.geometry,nothing,nothing)
 	V = A.fop(args)
-	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
+	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiWavefield):",A.name,typeof(A),eltype(V)]," / "))
 	return V
 end
 
@@ -161,7 +161,7 @@ function *(A::judiPDEadjoint{ADDT,ARDT},v::judiWavefield{vDT}) where {ADDT,ARDT,
 	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiWavefield):",A.name,typeof(A),vDT]," / "))
 	args = (A.geometry,nothing,nothing,v.data,nothing)
 	V = A.fop(args)
-	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
+	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiWavefield):",A.name,typeof(A),eltype(V)]," / "))
 	return V
 end