Emlink refactor (#2)

* major update. numeric completely working, priors added for expectation maximization

* removing temp files

* addendum. removing unnecessary temp file from emacs... again

.gitignore

@@ -23,5 +23,7 @@ docs/site/
 # environment.
 Manifest.toml
 
-
+actual_scratch.jl
+scratch.jl
+.#*
 \#*

Project.toml

@@ -1,19 +1,21 @@
 name = "FastLink"
 uuid = "11f39cfd-5548-489f-be9a-f4ad0ff6eadc"
 authors = ["Jack R. Williams <contact@jackryanwilliams.com>"]
-version = "0.1.1"
+version = "0.0.2"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 PooledArrays = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
 StringDistances = "88034a9c-02f8-509d-84a9-84ec65e18404"
 
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
 
 [targets]
 test = ["Test", "CSV", "Pkg"]

README.md

@@ -1,4 +1,7 @@
 # FastLink.jl
 Fast Probabilistic Record Linkage for the Julia Language
+## What is FastLink.jl
+
+The purpose of FastLink.jl is to bring a fast record linkage package to the julia language. When attempting to match large datasets using existing libraries in R and Python, I found they can be very slow and succumb to issues with memory pressure. This implementation of the fastlink algorithm is intended to scale effeciently in parallel and be able to easily handle matches between tabular data that span millions of rows. 
 
 [![Run tests](https://github.com/jw2249a/FastLink.jl/actions/workflows/test.yml/badge.svg)](https://github.com/jw2249a/FastLink.jl/actions/workflows/test.yml)

scratch.jl

@@ -1,108 +1,48 @@
 using Pkg
-Pkg.develop(path=".")
+#Pkg.develop(path=".")
 using DataFrames
 using BenchmarkTools
 using CSV
 using FastLink
 using PooledArrays
+import Pkg.Artifacts: @artifact_str
 
-numeric=false
-# files for performance
-test=true
-if test
-    a_fil="../dfA.csv"
-    b_fil="../dfB.csv"
-    if numeric
-        varnames=["housenum"]
-        match_method=["float"]
-        cut_a=[1]
-        cut_p=[2]
-    else
-        varnames=["firstname","middlename", "lastname","housenum"]
-        match_method=["string", "string","string", "float"]
-        cut_a=[0.92,0.92,0.92,1]
-        cut_p=[0.88,0.88,0.88,2]
-    end
-else
-    a_fil="../../rstudio/test_merge/data/test_a.csv"
-    b_fil="../../rstudio/test_merge/data/test_b.csv"
 
-    if numeric
-        varnames=["ZIP", "DOB_YEAR", "ZIP4"]
-        match_method=["float", "float", "float"]
-        cut_a=[1,1,1]
-        cut_p=[2,2,2]
-    else
-        varnames=["FIRST_NAME", "MIDDLE_NAME", "LAST_NAME", "STREET_NAME"]
-        cut_a=[0.92,0.92,0.92,0.92]
-        cut_p=[0.88,0.88,0.88,0.88]
-        #varnames=["FIRST_NAME", "MIDDLE_NAME", "LAST_NAME", "STREET_NAME", "STATE"]
-    end
-end
+a_fil = @artifact_str "dfA"
+b_fil = @artifact_str "dfB"
 
+varnames=["firstname","middlename", "lastname","housenum"]
+match_method=["string", "string","string", "float"]
+cut_a=[0.92,0.92,0.92,1]
+cut_p=[0.88,0.88,0.88,2]
 
-#[100,200,500,1_000,2_000,4_000, 5_000, 10_000,20_000, 40_000, 50_000,100_000,1_000_000]
-N1=10_000
-N2=500_000
 
 
-if test
-    dfA=CSV.read(a_fil, DataFrame,
-                 ntasks=1,
-                 pool=true,
-                 missingstring=["", "NA"])
-    dfB=CSV.read(b_fil, DataFrame,
-                 ntasks=1,
-                 pool=true,
-                 missingstring=["", "NA"])
-else
-    dfA=CSV.read(a_fil, DataFrame,
-                 limit=N1,
-                 ignoreemptyrows=true,
-                 ntasks=1,
-                 pool=true,
-                 missingstring=["", "NA", "NaN", "NULL", "Null"])
-    dfB=CSV.read(b_fil, DataFrame,
-                 limit=N2,
-                 ignoreemptyrows=true,
-                 ntasks=1,
-                 pool=true,
-                 missingstring=["", "NA", "NaN", "NULL", "Null"])
-end
+dfA=CSV.read("$(a_fil)/dfA.csv", DataFrame,
+             ntasks=1,
+             pool=true,
+             missingstring=["", "NA"])
+dfB=CSV.read("$(b_fil)/dfB.csv", DataFrame,
+             ntasks=1,
+             pool=true,
+             missingstring=["", "NA"])
 
 
-if !test && numeric
-    for var in varnames
-        dfA[!,var]=passmissing(x-> try return parse(Float64,x) catch e return 0.0 end).(dfA[:,var])
-        dfB[!,var]=passmissing(x-> try return parse(Float64,x) catch e return 0.0 end).(dfB[:,var])
-    end
+for var in varnames[1:3]
+    dfA[!,var] = PooledArray(passmissing(x->uppercase(x)).(dfA[:,var]))
+    dfB[!,var] = PooledArray(passmissing(x->uppercase(x)).(dfB[:,var]))
 end
 
-# if test && !numeric
-#     for var in varnames
-#         dfA[!,var] = PooledArray(passmissing(x->uppercase(x)).(dfA[:,var]))
-#         dfB[!,var] = PooledArray(passmissing(x->uppercase(x)).(dfB[:,var]))
-#     end
-# end
-
-
 config = fastLink(dfA,dfB,varnames,match_method=match_method,cut_a=cut_a,cut_p=cut_p,
                  threshold_match = 0.85)
 
-
-
 dump(config.fastlink_settings.comparison_funs[4])
 
 results=fastLink(dfA,dfB,varnames,match_method=match_method,cut_a=cut_a,cut_p=cut_p,
-
-
-
                  threshold_match = 0.85)()
 
-
-
-
 x=results[1].patterns_w
 x[findall(ismissing.(x.gamma_4) .== false .&& x.gamma_4 .== 1),:]
 x[findall(ismissing.(x.gamma_4)),:]
-44+7+1+43+79+1
+
+

src/DiBitMatrix.jl

@@ -0,0 +1,75 @@
+module DiBitMat
+import Base: getindex, setindex!, view
+import DataStructures: DiBitVector
+export DiBitMatrix
+
+"""
+Extending DiBitVectors from DataStructures.jl to include matrices.
+"""
+struct DiBitMatrix
+    data::DiBitVector
+    nrows::Integer
+    ncols::Integer
+end
+
+# base definition of the DiBitMatrix
+function DiBitMatrix(nrows::Integer, ncols::Integer)
+    data = DiBitVector(nrows * ncols, 0)  # Or choose an appropriate type
+    return DiBitMatrix(data, nrows, ncols)
+end
+
+# getting items by index
+function getindex(vm::DiBitMatrix, i::Int, j::Int)
+    linear_index = (j - 1) * vm.nrows + i
+    return vm.data[linear_index]
+end
+
+function getindex(vm::DiBitMatrix, ::Colon, j::Int)
+    column = zeros(UInt8, vm.nrows)
+    for i in 1:vm.nrows
+        linear_index = (j - 1) * vm.nrows + i
+        column[i] = vm.data[linear_index]
+    end
+    return column
+end
+
+function getindex(vm::DiBitMatrix, i::Int, ::Colon)
+    row = zeros(UInt8, vm.ncols)
+    for j in 1:vm.ncols
+        linear_index = (j - 1) * vm.nrows + i
+        row[j] = vm.data[linear_index]
+    end
+    return row
+end
+
+# setting items by index
+function setindex!(vm::DiBitMatrix, value::UInt8, i::T, j::T) where {T<:Integer}
+    linear_index = (j - 1) * vm.nrows + i
+    vm.data[linear_index] = value
+end
+
+# extending view to handle DiBitMatrix columns
+function getIndices(vm::DiBitMatrix,::Colon,j::Int)
+    return (j - 1) * vm.nrows + 1, (j - 1) * vm.nrows + vm.nrows
+end
+
+function getIndices(vm::DiBitMatrix, i::Int, ::Colon)
+    row = zeros(Integer, vm.ncols)
+    for j in 1:vm.ncols
+        row[j] = (j - 1) * vm.nrows + i
+    end
+    return row
+end
+
+function view(vm::DiBitMatrix,::Colon, j::Int)
+    start,finish=getIndices(vm,:,j)
+    return view(vm.data, start:finish)
+end
+
+function view(vm::DiBitMatrix, i::Int,::Colon)
+    vals=getIndices(vm, i,:)
+    return view(vm.data, vals)
+end
+
+
+end

src/FastLink.jl

@@ -3,17 +3,22 @@ using DataFrames
 import PooledArrays: PooledVector
 import Distributions: Dirichlet,rand
 
+# match constants
+const nonmatch::UInt8 = UInt8(0)
+const match1::UInt8 = UInt8(1)
+const match2::UInt8 = UInt8(2)
+const missingval::UInt8 = UInt8(3)
 
-include("resultMatrix.jl")
+include("DiBitMatrix.jl")
+using .DiBitMat
 include("gammas/Gammas.jl")
 using .Gammas
 
-include("tableCounts.jl")
+include("matchPatterns.jl")
 include("emlink.jl")
 include("getMatches.jl")
 include("fastlink/fastlink.jl")
 
-export(tableCounts)
 export(fastLink)
 
 

src/emlink.jl

@@ -24,33 +24,62 @@ end
 """
 Expectation maximization function. 
 """
-function emlinkMARmov(patterns::Dict, obs_a::Int,obs_b::Int,varnames::Vector{String}, ranges::Vector{UnitRange{Int64}}; p_m=0.1,iter_max=5000,tol=Float64(1e-05),missingval = [false,true])
+function emlinkMARmov(patterns::MatchPatterns, obs_a::Int, obs_b::Int,varnames::Vector{String};
+                      p_m=0.1,iter_max=5000,tol=1e-05, prior_lambda=0.0, w_lambda=0.0,
+                      prior_pi=0.0,w_pi=0.0, address_field=Vector{Bool}())
     # Initialize count and delta for while loop and break point
     delta = Float64(1)
     count = 1
 
     # Info for EM algorithm
     p_u = 1 - p_m
     nfeatures=length(varnames)
-    gamma_jk=collect(keys(patterns))
-    n_j = collect(values(patterns))
+    gamma_jk=patterns.patterns
+    n_j = length.(patterns.indices)
     N = length(n_j)
 
-    # TODO: add "if statement" λ priors are declared
-    psi = 1
-    mu = 1
-
-    ###########################################
-    # # TODO: add "if statement" for π priors #
-    # ## for address                          #
-    # ⍺₀_address = 1                          #
-    # ⍺₁_address = 1                          #
-    # address_field = falses(nfeatures)       #
-    # ## for lambda                           #
-    # ⍺₀_gender = 1                           #
-    # ⍺₁_gender = 1                           #
-    # genderaddress_field = falses(nfeatures) #
-    ###########################################
+    # if λ priors are declared
+    if prior_lambda == 0
+        psi = 1
+        mu = 1        
+    else
+        if w_lambda == 0
+            @error "If declaring a lambda prior, you need to declare weights via w_lambda."
+        elseif w_lambda > 0 | w_lambda < 0
+            @error "w_lambda must be between 0 and 1."
+        elseif w_lambda == 1
+            w_lambda = 1 - 1e-05
+        end
+        c_lambda = w_lambda/(1-w_lambda)
+        # hyperparameters for lambda
+        mu = prior_lambda * c_lambda * obs_a * obs_b + 1
+        psi = (1 - prior_lambda) * mu / prior_lambda
+    end
+
+    # if pi prior is declared
+    if prior_pi == 0
+        alpha0_address = 1
+        alpha1_address = 1
+        address_field = falses(nfeatures)
+    else
+        if prior_lambda == 0
+            @error "If declaring a prior on pi, you need to declare lambda prior."
+        elseif w_pi == 0
+            @error "If providing a prior for pi, please specify the weight using w_pi"
+        elseif w_pi < 0 | w_pi > 1
+            @error "w_pi must be between 0 and 1."
+        elseif w_pi == 1
+            w_pi = 1 - 1e-05
+        end
+
+        c_pi = w_pi / (1 - w_pi)
+        exp_match = prior_lambda * obs_a * obs_b
+
+        # Optimal hyperparameters for pi
+        alpha0_address = c_pi * prior_pi * exp_match + 1
+        alpha1_address = alpha0_address * (1 - prior_pi) / prior_pi
+    end
+
     # initialize variables that need value to be returned
     zeta_j=0.0
     num_prod = zeros(Float64,0)
@@ -65,8 +94,7 @@ function emlinkMARmov(patterns::Dict, obs_a::Int,obs_b::Int,varnames::Vector{Str
     p_gamma_kjm = missings(Union{Missing,Float64}, (nfeatures,N))
     p_gamma_kju = missings(Union{Missing,Float64}, (nfeatures,N))
     for c in 1:nfeatures
-        col=ranges[c]
-        vals_gamma_jk[c] = [i[col] == missingval ? missing : sum(i[col]) for i in gamma_jk]
+        vals_gamma_jk[c] = [i[c] == missingval ? missing : sum(i[c]) for i in gamma_jk]
         uvals_gamma_jk[c] = sort(unique([i for i in vals_gamma_jk[c] if !ismissing(i)]))
         c_m = collect(1:50:(length(uvals_gamma_jk[c])*50))
         p_gamma_km[c] = sort(rand(Dirichlet(c_m),1)[:],rev=false)
@@ -89,9 +117,13 @@ function emlinkMARmov(patterns::Dict, obs_a::Int,obs_b::Int,varnames::Vector{Str
         p_u = 1-p_m
 
         for i in 1:nfeatures
+            km_prob=sort([sum(num_prod[findall(skipmissing_equality(vals_gamma_jk[i], uvals_gamma_jk[i][j]))])
+                      for j in 1:length(uvals_gamma_jk[i])],rev=false)
+            if address_field[i]
+                km_prob += append!([alpha0_address], [alpha1_address for i in 1:(length(uvals_gamma_jk[i])-1)])
+            end
             p_gamma_km[i] = 
-                sort(probability_vector([sum(num_prod[findall(skipmissing_equality(vals_gamma_jk[i], uvals_gamma_jk[i][j]))])
-                      for j in 1:length(uvals_gamma_jk[i])]),rev=false)
+                probability_vector(km_prob)
             p_gamma_ku[i] =
                 sort(probability_vector([let sub1 = sub=findall(skipmissing_equality(vals_gamma_jk[i], uvals_gamma_jk[i][j]));
                                              sum(n_j[sub] - num_prod[sub])
@@ -118,6 +150,7 @@ function emlinkMARmov(patterns::Dict, obs_a::Int,obs_b::Int,varnames::Vector{Str
             pgamma_jm = p_gamma_jm, pgamma_ju = p_gamma_ju,
             patterns_w = data_w,
             patterns_b = gamma_jk,
+            indices = patterns.indices,
             iter_converge = count,
             obs_a = obs_a, obs_b = obs_b,
             varnames = varnames)