Various 20260711 : streamline read methods, reduce allocations

.github/workflows/ci.yml

@@ -15,8 +15,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
-          - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
+          - '1.9' # minimum Julia version supported
+          - '1.10' # current LTS version of Julia
+          - '1' # latest stable 1.x release of Julia.
 #          - 'nightly'
         os:
           - ubuntu-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "MeshArrays"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
-version = "0.5.7"
+version = "0.5.8"
 authors = ["gaelforget <gforget@mit.edu>"]
 
 [deps]

src/GridPaths.jl

@@ -47,7 +47,7 @@ function rotate_XCYC(Γ,R)
 
 	tmpx2=tmp3[1,:]; tmpy2=tmp3[2,:]; tmpz2=tmp3[3,:]
 	tmpx[tmp1]=tmpx2; tmpy[tmp1]=tmpy2; tmpz[tmp1]=tmpz2
-	x=γ.read(tmpx,Γ.XC); y=γ.read(tmpy,Γ.XC); z=γ.read(tmpz,Γ.XC)
+	read!(tmpx,x); read!(tmpy,y); read!(tmpz,z)
 
 	x,y,z
 end
@@ -75,15 +75,15 @@ function shorter_paths!(xyz,xyz0,msk_in)
 	msk_out=[]
 	for kk in 1:3
 		#select field to treat:
-		mm=msk_in[kk]
+		mm=similar(msk_in[kk])
 		#select the shorther segment:
-		tmpm=γ.write(mm)
+		tmpm=γ.write(msk_in[kk])
        if theta[2]-theta[1]<=pi
             tmpm[findall( (tmptheta.>theta[2]).|(tmptheta.<theta[1]) )].=0.0
         else
             tmpm[findall( (tmptheta.<=theta[2]).&(tmptheta.>=theta[1]) )].=0.0
         end
-        mm=γ.read(tmpm,mm)
+        read!(tmpm,mm)
         #store result:
 		push!(msk_out,mm)
     end

src/ReadWrite.jl

@@ -13,27 +13,37 @@ function nFacesEtc(a::AbstractMeshArray)
   return nFaces, n3, n4
 end
 
+function _load_binary(fil::String, y::AbstractMeshArray)
+  (n1,n2) = y.grid.ioSize
+  (_,n3,n4) = nFacesEtc(y)
+  fid = open(fil)
+  xx = Array{eltype(y),1}(undef, n1*n2*n3*n4)
+  read!(fid, xx)
+  xx = reshape(hton.(xx), (n1,n2,n3,n4))
+  close(fid)
+  xx
+end
+
 """
     read(fil::String,x::AbstractMeshArray)
 
-Read array from file and return as a MeshArray. 
+Read array from file and return as a MeshArray.
 
 _The second argument (MeshArray or gcmgrid) provides the grid specifications (x.grid.ioSize)._
 ```
 """
-function read(fil::String,x::AbstractMeshArray)
-
-  (n1,n2)=x.grid.ioSize
-  (nFaces,n3,n4)=nFacesEtc(x)
+function read(fil::String, x::AbstractMeshArray)
+  y = similar(x; m=x.meta)
+  read!(fil, y)
+  return y
+end
 
-  fid = open(fil)
-  xx = Array{eltype(x),1}(undef,n1*n2*n3*n4)
-  read!(fid,xx)
-  xx = reshape(hton.(xx),(n1,n2,n3,n4))
-  close(fid)
+"""
+    read!(fil::String, y::AbstractMeshArray)
 
-  return x.grid.read(xx,x)
-end
+Read binary file `fil` into pre-allocated MeshArray `y` (in-place, no allocation).
+"""
+read!(fil::String, y::AbstractMeshArray) = read!(_load_binary(fil, y), y)
 
 """
     read(xx::Array,γ::gcmgrid)
@@ -63,12 +73,14 @@ function read(xx::Array,γ::gcmgrid; verbose=false)
   verbose ? println(size(yy)) : nothing
 
   if n3==1&&n4==1
-    read(yy,MeshArray(γ,γ.ioPrec))
+    y = MeshArray(γ,γ.ioPrec)
   elseif n4==1
-    read(yy,MeshArray(γ,γ.ioPrec,n3))
+    y = MeshArray(γ,γ.ioPrec,n3)
   else
-    read(yy,MeshArray(γ,γ.ioPrec,n3,n4))
+    y = MeshArray(γ,γ.ioPrec,n3,n4)
   end
+  read!(yy,y)
+  return y
 end
 
 """
@@ -95,7 +107,7 @@ function read!(xx::Array,x::AbstractMeshArray)
   tmp=zeros(x.grid)
   for i3 in 1:n3
     for i4 in 1:n4
-      read_one!(xx[:,:,i3,i4],tmp)
+      _read_one!(xx[:,:,i3,i4],tmp)
       for f in 1:nFaces
         if (n3>1)&&(n4>1)
           x[f,i3,i4].=tmp[f]
@@ -110,11 +122,11 @@ function read!(xx::Array,x::AbstractMeshArray)
 end
 
 """
-    read!(xx::Array,x::AbstractMeshArray)
+    _read_one!(xx::Array,x::AbstractMeshArray)
 
-Reformat one array of size x.grid.ioSize, and write **in-place** into MeshArray `x``.
+Reformat one array of size x.grid.ioSize, and write **in-place** into MeshArray `x`.
 """
-function read_one!(xx::Array,x::AbstractMeshArray; verbose=false)
+function _read_one!(xx::Array,x::AbstractMeshArray; verbose=false)
   test1=in(x.grid.class,["PeriodicChannel","PeriodicDomain"])
   format=(test1 ? :simple : :compact)
 
@@ -128,12 +140,16 @@ function read_one!(xx::Array,x::AbstractMeshArray; verbose=false)
     if format==:compact
       i0=i1+1
       i1=i1+nn*mm
-      x.f[iFace]=reshape(xx[:][i0:i1,:],(nn,mm))
+      for jj in 1:mm, ii in 1:nn
+        x.f[iFace][ii,jj] = xx[i0 + (ii-1) + (jj-1)*nn]
+      end
     else
       i0=(mod(i1,n1)==0 ? 1 : i1+1)
       j0=(mod(i1,n1)==0&&iFace!==1 ? j0+mm : j0)
       i1=i0+nn-1; j1=j0+mm-1
-      x.f[iFace]=reshape(xx[i0:i1,j0:j1],(nn,mm))
+      for jj in 1:mm, ii in 1:nn
+        x.f[iFace][ii,jj] = xx[i0+ii-1, j0+jj-1]
+      end
     end
 end
 

src/grids/main.jl

@@ -210,7 +210,9 @@ function GridLoadVar(nam::String,γ::gcmgrid)
     if sum(nam.==list_n)==1
         ii=findall(nam.==list_n)[1]
         m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii])
-        tmp1=γ.read(joinpath(γ.path,list_n[ii]*".data"),MeshArray(γ,γ.ioPrec;meta=m))
+        tmp1=MeshArray(γ,γ.ioPrec;meta=m)
+        read!(joinpath(γ.path,list_n[ii]*".data"),tmp1)
+        return tmp1
     elseif sum(nam.==list1d_n)==1
         fil=joinpath(γ.path,nam*".data")
         γ.ioPrec==Float64 ? reclen=8 : reclen=4
@@ -220,16 +222,19 @@ function GridLoadVar(nam::String,γ::gcmgrid)
         tmp1 = Array{γ.ioPrec,1}(undef,n3)
         read!(fid,tmp1)
         tmp1 = hton.(tmp1)
+        return tmp1
     elseif sum(nam.==list3d_n)==1
         fil=joinpath(γ.path,"RC.data")
         γ.ioPrec==Float64 ? reclen=8 : reclen=4
         n3=Int64(stat(fil).size/reclen)
 
         ii=findall(nam.==list3d_n)[1]
         m=varmeta(list3d_u[ii],list3d_p[ii],missing,list3d_n[ii],list3d_n[ii])
-        tmp1=γ.read(joinpath(γ.path,list3d_n[ii]*".data"),MeshArray(γ,γ.ioPrec,n3;meta=m))
+        tmp1=MeshArray(γ,γ.ioPrec,n3;meta=m)
+        read!(joinpath(γ.path,list3d_n[ii]*".data"),tmp1)
+        return tmp1
     else
-        tmp1=missing
+        return missing
     end
 end
 

src/grids/simple.jl

@@ -67,7 +67,7 @@ function GridLoad_ones(γ::gcmgrid;option="minimal")
 
     pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.]
     if option=="full"
-        list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","hFacC","hFacS","hFacW")
+        list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","hFacC","hFacW","hFacS")
         list_u=(u"m",u"m",u"m",u"m",u"m^2",u"m^2",u"m^2",u"m^2",u"m",u"m",u"m",u"m",u"m",1.0,1.0,1.0)
         list_p=(pc,pg,pc,pg,pc,pu,pv,pg,pu,pv,pv,pu,pc,fill(0.5,3),[0.,0.5,0.5],[0.5,0.,0.5])
     elseif option=="light"
@@ -79,12 +79,13 @@ function GridLoad_ones(γ::gcmgrid;option="minimal")
         list_u=(u"°",u"°")
         list_p=(pc,pc)
     end
-
+
+    tmp1=fill(1.,(ioSize[:]))
     for ii=1:length(list_n)
-        tmp1=fill(1.,(ioSize[:]))
         m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii])
-        tmp1=γ.read(tmp1,MeshArray(γ,Float64;meta=m))
-        Γ[list_n[ii]]=tmp1
+        tmp2=MeshArray(γ,Float64;meta=m)
+        read!(tmp1,tmp2)
+        Γ[list_n[ii]]=tmp2
     end
 
     XC=[i-0.5 for i in 1:ioSize[1], j in 1:ioSize[2]]
@@ -94,7 +95,7 @@ function GridLoad_ones(γ::gcmgrid;option="minimal")
         YG=[j-1.0 for i in 1:ioSize[1], j in 1:ioSize[2]]
     end
 
-    (mp,mq)=Int.(γ.ioSize[:]./γ.fSize[1][:])
+    (mp,mq)=Int.(γ.ioSize[:]./γ.fSize[1][:]) #wont work with LLC grid
     for i in 1:nFaces
         (np,nq)=γ.fSize[i]
         ip=(1:np) .+ (mod1(i,mp)-1).*np
@@ -104,7 +105,7 @@ function GridLoad_ones(γ::gcmgrid;option="minimal")
         Γ["YC"][i]=YC[ip,iq]
         option=="full" ? Γ["YG"][i]=YG[ip,iq] : nothing
     end
-    
+
     Dict_to_NamedTuple(Γ)
 end
 
@@ -121,17 +122,8 @@ etc accordingly with global `XC,YC` etc .
 function UnitGrid(ioSize::NTuple{2, Int},tileSize::NTuple{2, Int}; option="minimal")
     nF=div(prod(ioSize),prod(tileSize))
     fSize=fill(tileSize,nF)
-
     γ=gcmgrid("","PeriodicDomain",nF,fSize, ioSize, Float32, read_tiles, write_tiles)
     Γ=GridLoad_ones(γ;option=option)
-
-    Γ.XC[:]=γ.read([i-0.5 for i in 1:ioSize[1], j in 1:ioSize[2]],γ)
-    Γ.YC[:]=γ.read([j-0.5 for i in 1:ioSize[1], j in 1:ioSize[2]],γ)
-    if option=="full"
-        Γ.XG[:]=γ.read([i-1.0 for i in 1:ioSize[1], j in 1:ioSize[2]],γ)
-        Γ.YG[:]=γ.read([j-1.0 for i in 1:ioSize[1], j in 1:ioSize[2]],γ)
-    end
-
     return Γ,γ
 end
 

src/grids/tiles.jl

@@ -43,8 +43,12 @@ Return an `Array` of tiles which cover `x` according to tile partition `τ`.
 
 ```jldoctest; output = false
 using MeshArrays
-γ=GridSpec("LatLonCap",MeshArrays.Dataset("GRID_LLC90"))
-d=γ.read(γ.path*"Depth.data",MeshArray(γ,γ.ioPrec))
+
+γ=GridSpec(ID=:LLC90)
+d=MeshArray(γ,γ.ioPrec)
+f=joinpath(γ.path,"Depth.data")
+read!(f,d)
+
 τ=Tiles(γ,30,30)
 td=Tiles(τ,d)
 

src/types/gcmarray.jl

@@ -40,10 +40,11 @@ function gcmarray(A::Array{T,N};
   ios=[s[1] s[2]]
   γ=gcmgrid("","PeriodicDomain",1, fs, ios, T, read, write)
   if N==2
-    B=γ.read(A,gcmarray(γ,T;meta=meta))
+    B=gcmarray(γ,T;meta=meta)
   else
-    B=γ.read(A,gcmarray(γ,T,size(A)[3];meta=meta))
+    B=gcmarray(γ,T,size(A)[3];meta=meta)
   end
+  read!(A,B)
   B
 end
 

test/testsets/transport.jl

@@ -6,8 +6,8 @@
     @suppress show(Γ.XC)
 
     path=MeshArrays.Dataset("GRID_LLC90")
-    Tx=γ.read(joinpath(path,"TrspX.bin"),MeshArray(γ,Float32))
-    Ty=γ.read(joinpath(path,"TrspY.bin"),MeshArray(γ,Float32))
+    Tx=MeshArray(γ,Float32); read!(joinpath(path,"TrspX.bin"),Tx)
+    Ty=MeshArray(γ,Float32); read!(joinpath(path,"TrspY.bin"),Ty)
     plot(Γ.XC)
 
     hFacC=GridLoadVar("hFacC",γ)
@@ -56,8 +56,13 @@
     ones(y)
     zeros(y)
 
-    GM_PsiX=read(randn(90,1170,50),Γ.hFacW)
-    GM_PsiY=read(randn(90,1170,50),Γ.hFacS)
+    m=varmeta(1.0, [0.0, 0.5, 0.5], missing, "GM_PsiX", "GM_PsiX")
+    GM_PsiX=MeshArray(γ,γ.ioPrec,50, meta=m)
+    read!(randn(γ.ioSize...,50),GM_PsiX)
+    m=varmeta(1.0, [0.5, 0.0, 0.5], missing, "GM_PsiY", "GM_PsiY")
+    GM_PsiY=MeshArray(γ,γ.ioPrec,50, meta=m)
+    read!(randn(γ.ioSize...,50),GM_PsiY)
+
     bolusU, bolusV, bolusW=MeshArrays.calc_bolus(GM_PsiX,GM_PsiY, Γ)
 
     read(rand(90*1170),γ)