Merge pull request #275 from slimgroup/pythoncall

Switch to  Pythoncall

.github/workflows/ci-examples.yml

@@ -42,6 +42,7 @@ jobs:
       DEVITO_LANGUAGE: "openmp"
       DEVITO_LOGGING: "ERROR"
       OMP_NUM_THREADS: 1
+      JULIA_NUM_THREADS: 1
       NITER: 2
 
     strategy:
@@ -60,19 +61,23 @@ jobs:
           version: ${{ matrix.version }}
           arch: x64
 
+      - uses: julia-actions/cache@v2
+
       - name: Set julia python
         run: |
           python3 -m pip install -U pip
           python3 -m pip install "matplotlib<3.9" seiscm colorcet
+          echo "JULIA_PYTHONCALL_EXE=$(which python3)" >> $GITHUB_ENV
+          echo "JULIA_CONDAPKG_BACKEND=\"Null\"" >> $GITHUB_ENV
           echo "PYTHON=$(which python3)" >> $GITHUB_ENV
-          julia -e 'using Pkg;Pkg.add(["PyCall", "PyPlot", "SlimPlotting"])'
+          echo "PYCALL_JL_RUNTIME_PYTHON=$(which python3)" >> $GITHUB_ENV
 
       - name: Build JUDI
         uses: julia-actions/julia-buildpkg@latest
 
       - name: Install packages
         run: |
-          julia -e 'using Pkg;Pkg.add(["NLopt", "Flux", "JOLI", "Zygote", "IterativeSolvers", "SlimOptim", "HDF5", "SegyIO", "SetIntersectionProjection"])'
+          julia -e 'using Pkg;Pkg.add(["SlimPlotting", "PyPlot", "NLopt", "Flux", "JOLI", "Zygote", "IterativeSolvers", "SlimOptim", "HDF5", "SegyIO", "SetIntersectionProjection"])'
           julia -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd()))'
 
       - name: Run examples

.github/workflows/ci-judi.yml

@@ -22,6 +22,7 @@ jobs:
       DEVITO_ARCH: ${{ matrix.cc }}
       DEVITO_LANGUAGE: "openmp"
       OMP_NUM_THREADS: 4
+      JULIA_NUM_THREADS: 1
       GROUP: "JUDI"
 
     strategy:
@@ -49,23 +50,13 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
 
+      - uses: julia-actions/cache@v2
+
       - name: Setup clang for osx
         if: runner.os == 'macOS'
         run: |
           brew install llvm libomp
-          echo "/opt/homebrew/bin:/opt/homebrew/opt/llvm/bin" >> $GITHUB_PATH
-
-      - name: Set up Python 3.9
-        uses: actions/setup-python@v5
-        with:
-          python-version: 3.9
-
-      - name: Set julia python
-        run: |
-          echo "PYTHON=$(which python3)" >> $GITHUB_ENV
-          echo "PYCALL_JL_RUNTIME_PYTHON=$(which python3)" >> $GITHUB_ENV
-          python3 -m pip install devito[tests,extras]@git+https://github.com/devitocodes/devito.git
-          PYCALL_JL_RUNTIME_PYTHON=$(which python3) PYTHON=$(which python3) julia -e 'using Pkg;Pkg.add("PyCall");Pkg.build("PyCall")'
+          echo "/opt/homebrew/opt/llvm/bin" >> $GITHUB_PATH
 
       - name: Build JUDI
         uses: julia-actions/julia-buildpkg@latest

.github/workflows/ci-op.yml

@@ -23,6 +23,7 @@ jobs:
       DEVITO_LANGUAGE: "openmp"
       DEVITO_LOGGING: "INFO"
       OMP_NUM_THREADS: ${{ matrix.omp }}
+      JULIA_NUM_THREADS: 1
       GROUP: ${{ matrix.op }}
 
     strategy:
@@ -82,11 +83,13 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
 
+      - uses: julia-actions/cache@v2
+
       - name: Setup clang for osx
         if: runner.os == 'macOS'
         run: |
           brew install llvm libomp
-          echo "/opt/homebrew/bin:/opt/homebrew/opt/llvm/bin" >> $GITHUB_PATH
+          echo "/opt/homebrew/opt/llvm/bin" >> $GITHUB_PATH
 
       - name: Set up Python 3.9
         uses: actions/setup-python@v5
@@ -97,8 +100,8 @@ jobs:
         run: |
           echo "PYTHON=$(which python3)" >> $GITHUB_ENV
           echo "PYCALL_JL_RUNTIME_PYTHON=$(which python3)" >> $GITHUB_ENV
-          python3 -m pip install devito[tests,extras]@git+https://github.com/devitocodes/devito.git
-          PYCALL_JL_RUNTIME_PYTHON=$(which python3) PYTHON=$(which python3) julia -e 'using Pkg;Pkg.add("PyCall");Pkg.build("PyCall")'
+          echo "JULIA_PYTHONCALL_EXE=$(which python3)" >> $GITHUB_ENV
+          echo "JULIA_CONDAPKG_BACKEND=\"Null\"" >> $GITHUB_ENV
 
       - name: Build JUDI
         uses: julia-actions/julia-buildpkg@latest

.github/workflows/deploy_doc.yaml

@@ -30,8 +30,11 @@ jobs:
 
       - name: Set julia python
         run: |
-          PYTHON=$(which python3) julia -e 'using Pkg;Pkg.add("PyCall");Pkg.build("PyCall")'
-  
+          echo "PYTHON=$(which python3)" >> $GITHUB_ENV
+          echo "PYCALL_JL_RUNTIME_PYTHON=$(which python3)" >> $GITHUB_ENV
+          echo "JULIA_PYTHONCALL_EXE=$(which python3)" >> $GITHUB_ENV
+          echo "JULIA_CONDAPKG_BACKEND=\"Null\"" >> $GITHUB_ENV
+
       - name: Build JUDI
         uses: julia-actions/julia-buildpkg@latest
 

Project.toml

@@ -1,7 +1,7 @@
 name = "JUDI"
 uuid = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
 authors = ["Philipp Witte, Mathias Louboutin"]
-version = "3.4.7"
+version = "4.0.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -13,7 +13,7 @@ JOLI = "bb331ad6-a1cf-11e9-23da-9bcb53c69f6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SegyIO = "157a0f19-4d44-4de5-a0d0-07e2f0ac4dfa"
@@ -25,9 +25,9 @@ DSP = "0.6, 0.7"
 Dierckx = "0.4, 0.5"
 FFTW = "1"
 Flux = "0.12, 0.13, 0.14"
-JOLI = "0.7, 0.8"
+JOLI = "0.9"
 OrderedCollections = "1.5"
-PyCall = "1.18, 1.90, 1.91, 1.62"
+PythonCall = "0.9"
 Requires = "1"
 SegyIO = "0.7.7 - 0.8.5"
 TimerOutputs = "0.5"
@@ -54,4 +54,4 @@ test = ["Aqua", "JLD2", "Printf", "Test", "TimerOutputs", "Flux"]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

deps/build.jl

@@ -1,21 +1,26 @@
-using PyCall
+using PythonCall
 
 # Check devito version and update if necessary
 struct DevitoException <: Exception
     msg::String
 end
 
+pyexe = PythonCall.python_executable_path()
+
 pk = try
     pyimport("pkg_resources")
 catch e
-    run(PyCall.python_cmd(`-m pip install -U --user --no-cache-dir setuptools`))
+    run(Cmd(`$(pyexe) -m pip install -U --user --no-cache-dir setuptools`))
     pyimport("pkg_resources")
 end
 
+################## JOLI ##################
+run(Cmd(`$(pyexe) -m pip install -U --user --no-cache-dir PyWavelets`))
+
 ################## Devito ##################
 # pip command
 dvver = "4.8.10"
-cmd = PyCall.python_cmd(`-m pip install --user --no-cache-dir devito\[extras,tests\]\>\=$(dvver)`)
+cmd = Cmd(`$(pyexe) -m pip install --user --no-cache-dir devito\[extras,tests\]\>\=$(dvver)`)
 
 try
     dv_ver = VersionNumber(split(pk.get_distribution("devito").version, "+")[1])
@@ -24,13 +29,13 @@ try
         run(cmd)
     end
 catch e
-    @info "Devito  not installed, installing with PyCall python"
+    @info "Devito  not installed, installing with PythonCall python"
     run(cmd)
 end
 
 ################## Matplotlib ##################
 try
     mpl = pyimport("matplotlib")
 catch e
-    run(PyCall.python_cmd(`-m pip install --user --no-cache-dir matplotlib`))
+    run(Cmd(`$(pyexe) -m pip install --user --no-cache-dir matplotlib`))
 end

examples/scripts/modeling_basic_2D.jl

@@ -42,7 +42,7 @@ m0 = (1f0 ./ v0).^2
 dm = vec(m0 - m)
 
 # Setup model structure
-nsrc = 2	# number of sources
+nsrc = 3	# number of sources
 model = Model(n, d, o, m)
 model0 = Model(n, d, o, m0)
 
@@ -121,7 +121,7 @@ dobs = Pr*F*adjoint(Ps)*q
 
 #' Plot the shot record
 fig = figure()
-plot_sdata(dobs[1]; new_fig=false, name="Synthetic data", cmap=dcmap)
+plot_sdata(dobs[2]; new_fig=false, name="Synthetic data", cmap=dcmap)
 display(fig)
 
 #' Because we have abstracted the linear algebra, we can solve the adjoint wave-equation as well 
@@ -137,7 +137,7 @@ dot1 = dot(q, qad)
 # <F x, y>
 dot2 = dot(dobs, dobs)
 # Compare
-@show dot1, dot2, (dot2 - dot2)/(dot1 + dot2)
+@show dot1, dot2, (dot1 - dot2)/(dot1 + dot2)
 
 #' # Inversion
 #' Our main goal is to provide an inversion framework for seismic inversion. To this end, as shown earlier,
@@ -151,15 +151,25 @@ rtm = adjoint(J)*dD
 
 #' We show the linearized data.
 fig = figure()
-plot_sdata(dobs[1]; new_fig=false, name="Linearized data", cmap=dcmap)
+plot_sdata(dD[2]; new_fig=false, name="Linearized data", cmap=dcmap)
 display(fig)
 
-
 #' And the RTM image
 fig = figure()
 plot_simage(rtm'; new_fig=false, name="RTM image", cmap=imcmap)
 display(fig)
 
+#' We can easily now again test the adjointness of our operator with the standard dot test. Because we
+#' intend to conserve our linear algebra abstraction, `judiVector` implements all the necessary linear 
+#' algebra functions such as dot product or norm to be used directly.
+# <x, J'y>
+dot3 = dot(dm, rtm)
+# <J x, y>
+dot4 = dot(dD, dD)
+# Compare
+@show dot3, dot4, (dot3 - dot4)/(dot3 + dot4)
+
+
 #' ## Inversion utility functions
 #' We currently introduced the lineaar operators that allow to write seismic modeling and inversion in a high-level, linear algebra way. These linear operators allow the script to closely follow the mathematics and to be readable and understandable.
 #' 
@@ -213,13 +223,16 @@ display(fig)
 #' Finally, JUDI implements TWRI, an augmented method to tackle cycle skipping. Once again we provide a computationnally efficient wrapper function that returns the objective value and necessary gradients
 f, gm, gy = twri_objective(model0, q, dobs, nothing; options=opt, optionswri=TWRIOptions(params=:all))
 # With on-the-fly DFT, experimental
-f, gmf = twri_objective(model0, q, dobs, nothing; options=Options(frequencies=[[.009, .011], [.008, .012]]), optionswri=TWRIOptions(params=:m))
+freqs = [[.009, .011], [.008, .012], [0.007, 0.0010]]
+f, gmf = twri_objective(model0, q, dobs, nothing; options=Options(frequencies=freqs), optionswri=TWRIOptions(params=:m))
 
 #' Plot gradients
 fig = figure()
 plot_simage(gm'; new_fig=false, name="TWRI gradient w.r.t m", cmap=imcmap)
 display(fig)
 
 fig = figure()
-plot_sdata(gy[1]; new_fig=false, name="TWRI gradient w.r.t y", cmap=dcmap)
-display(fig)
+plot_sdata(gy[2]; new_fig=false, name="TWRI gradient w.r.t y", cmap=dcmap)
+display(fig)
+
+println("All done")

ext/JLD2JUDIExt.jl

@@ -10,7 +10,7 @@ JUDI.Geometry(x::JLD2.ReconstructedMutable{N, FN, NT}) where {N, FN, NT} = Geome
 function JUDI.tof32(x::JLD2.ReconstructedStatic{N, FN, NT}) where {N, FN, NT}
     # Drop "typed" signature
     reconstructT = Symbol(split(string(N), "{")[1])
-    return JUDI.tof32(eval(reconstructT)([getproperty(x, f) for f in FN]...))
+    return JUDI.tof32(getproperty(@__MODULE__, reconstructT)([getproperty(x, f) for f in FN]...))
 end
 
 end

src/JUDI.jl

@@ -18,7 +18,7 @@ end
 using LinearAlgebra, Random, Printf
 using Distributed
 using DSP, FFTW, Dierckx
-using PyCall
+using PythonCall
 using JOLI, SegyIO
 using ChainRulesCore
 using OrderedCollections
@@ -50,19 +50,17 @@ import JOLI: jo_convert
 # FFTW
 import FFTW: fft, ifft
 
-# Import pycall array to python for easy plotting
-import PyCall.NpyArray
-
 # Import AD rrule
 import ChainRulesCore: rrule
 
 # Set python paths
 export devito
 
-const pm = PyNULL()
-const ac = PyNULL()
-const pyut = PyNULL()
-const devito = PyNULL()
+const pm = PythonCall.pynew()
+const ac = PythonCall.pynew()
+const pyut = PythonCall.pynew()
+const devito = PythonCall.pynew()
+const np = PythonCall.pynew()
 
 # Constants
 nworkers(::Any) = length(workers())
@@ -99,13 +97,12 @@ include("rrules.jl")
 
 # Initialize
 function __init__()
-    pushfirst!(PyVector(pyimport("sys")."path"), joinpath(JUDIPATH, "pysource"))
-    copy!(pm, pyimport("models"))
-    copy!(ac, pyimport("interface"))
-    copy!(pyut, pyimport("utils"))
-    copy!(devito, pyimport("devito"))
-    # Initialize lock at session start
-    PYLOCK[] = ReentrantLock()
+    pyimport("sys").path.append(joinpath(JUDIPATH, "pysource"))
+    PythonCall.pycopy!(pm, pyimport("models"))
+    PythonCall.pycopy!(ac, pyimport("interface"))
+    PythonCall.pycopy!(pyut, pyimport("utils"))
+    PythonCall.pycopy!(devito, pyimport("devito"))
+    PythonCall.pycopy!(np, pyimport("numpy"))
 
     # Make sure there is no conflict for the cuda init thread with CUDA.jl
     if get(ENV, "DEVITO_PLATFORM", "") == "nvidiaX"

src/TimeModeling/LinearOperators/lazy.jl

@@ -215,4 +215,4 @@ _as_src(::judiNoopOperator, ::AbstractModel, q::judiMultiSourceVector) = q
 
 ############################################################################################################################
 ###### Evaluate lazy operation
-eval(rhs::judiRHS) = rhs.d
+eval_lazy(rhs::judiRHS) = rhs.d

src/TimeModeling/Modeling/losses.jl

@@ -12,14 +12,12 @@ and its derivative w.r.t `x`
     `x-y`
 
 """
-function mse(x::Matrix{T}, y::Matrix{T}) where {T<:Number}
-    f = .5f0 * norm(x - y, 2)^2
+function mse(x::AbstractArray{T}, y::AbstractArray{T}) where {T<:Number}
+    f = T(.5) * norm(x - y, 2)^2
     r = x - y
-    return f::T, r::Matrix{T}
+    return f, r
 end
 
-_studentst_loss(x::T, y::T, k::T) where {T<:Number} = T(1/2) * (k + 1) * log(1 + (x-y)^2 / k)
-
 """
 studentst(x, y)
 
@@ -32,9 +30,11 @@ and its derivative w.r.t x
     `(k + 1) * (x - y) / (k + (x - y)^2)`
 
 """
-function studentst(x::Matrix{T}, y::Matrix{T}; k=T(2)) where {T<:Number}
+function studentst(x::AbstractArray{T}, y::AbstractArray{T}; k=T(2)) where {T<:Number}
     k = convert(T, k)
     f = sum(_studentst_loss.(x, y, k))
     r = (k + 1) .* (x - y) ./ (k .+ (x - y).^2)
-    return f::T, r::Matrix{T}
-end
+    return f, r
+end
+
+_studentst_loss(x::T, y::T, k::T) where {T<:Number} = T(1/2) * (k + 1) * log(1 + (x-y)^2 / k)