Skip to content

Commit

Permalink
Initial
Browse files Browse the repository at this point in the history
  • Loading branch information
yufongpeng committed Nov 29, 2023

Verified

This commit was signed with the committer’s verified signature.
1 parent 5f0db23 commit bfa1d17
Showing 5 changed files with 317 additions and 1 deletion.
9 changes: 9 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -3,8 +3,17 @@ uuid = "3ed48883-6809-43ec-b77a-d23b0650d8c4"
authors = ["Yu-Fong Peng <sciphypar@gmail.com>"]
version = "1.0.0-DEV"

[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"

[compat]
julia = "1.9"
CSV = "0.10"
GLM = "1.9"
TypedTables = "1.4"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
58 changes: 57 additions & 1 deletion src/Calibration.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,61 @@
module Calibration

# Write your package code here.
using GLM, CSV, TypedTables, LinearAlgebra
export AbstractCalibration, MultipleCalibration, SingleCalibration,
AbstractAnalysisTable, ColumnAnalysisTable, RowAnalysisTable,
Project, project, calibration,
read_calibration, read_analysistable, read_project,
cal_range, lloq, uloq, accuracy, accuracy!,
inv_predict, inv_predict_cal!, inv_predict_accuracy!,
formula_repr, weight_repr, formula_repr_utf8, weight_repr_utf8, format_number

import Base: getproperty, show

abstract type AbstractCalibration{T} end
abstract type AbstractAnalysisTable{T} end

mutable struct MultipleCalibration{T} <: AbstractCalibration{T}
analyte::Int
isd::Int
type::Bool
zero::Bool
weight::Float64
formula::FormulaTerm
source::AbstractAnalysisTable{T}
table::Table
model
end

mutable struct SingleCalibration{T} <: AbstractCalibration{T}
analyte::Int
isd::Int
source::AbstractAnalysisTable{T}
β::Float64
end

struct ColumnAnalysisTable{T} <: AbstractAnalysisTable{T}
sample_name::Symbol
analyte_name::Vector{Symbol}
analytes::Vector{T}
table
end
getproperty(tbl::ColumnAnalysisTable, property::Symbol) = property == :samples ? Symbol.(getproperty(tbl.table, tbl.sample_name)) : getfield(tbl, property)

struct RowAnalysisTable{T} <: AbstractAnalysisTable{T}
sample_name::Vector{Symbol}
analyte_name::Symbol
analytes::Vector{T}
table
end
getproperty(tbl::RowAnalysisTable, property::Symbol) = property == :samples ? getfield(tbl, :sample_name) : getfield(tbl, property)

mutable struct Project
calibration::Vector{<: AbstractCalibration}
sample::Union{AbstractAnalysisTable, Nothing}
end

include("utils.jl")
include("cal.jl")
include("io.jl")

end
88 changes: 88 additions & 0 deletions src/cal.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,88 @@

inv_predict(project::Project, tbl) = inv_predict(project.calibration, tbl)
inv_predict_cal!(project::Project) = (inv_predict_cal!.(project.calibration); project)
inv_predict_cal!(cal::SingleCalibration) = cal
function inv_predict_cal!(cal::MultipleCalibration)
cal.table.x̂ .= inv_predict(cal, cal.table.y)
cal
end
inv_predict(cal::MultipleCalibration, tbl::AbstractAnalysisTable) = inv_predict(cal, get_analyte(tbl, cal.source.analytes[cal.analyte]))
inv_predict(cal::SingleCalibration, tbl::AbstractAnalysisTable) = inv_predict(cal, get_analyte(tbl, cal.source.analytes[cal.analyte]))
function inv_predict(cal::MultipleCalibration, y::AbstractArray)
β = cal.model.model.pp.beta0
if cal.type && cal.zero
y ./ β[1]
elseif cal.type
(y .- β[1]) ./ β[2]
else
c, b, a = cal.zero ? (0, β...) : β
d = @. max(b ^ 2 + 4 * a * (y - c), 0)
@. (-b + sqrt(d)) / 2a
end
end
inv_predict(cal::SingleCalibration, y::AbstractArray) = y ./ cal.β

accuracy(project::Project, tbl = project.calibration.table) = accuracy(project.calibration, tbl)
accuracy(cal::MultipleCalibration, tbl = cal.table.y) = accuracy(inv_predict(cal, tbl), tbl.x)
accuracy(cal::SingleCalibration, tbl) = accuracy(inv_predict(cal, tbl), tbl.x)
accuracy!(project::Project) = (accuracy!(project.calibration); project)
function accuracy!(cal::MultipleCalibration)
cal.table.accuracy .= accuracy(cal.table.x̂, cal.table.x)
cal
end
accuracy!(cal::SingleCalibration) = cal
accuracy(x̂::AbstractVector, x::AbstractVector) = @./ x

inv_predict_accuracy! = accuracy! inv_predict_cal!

function calibration(tbl::Table, source::AbstractAnalysisTable;
analyte = 1,
isd = 0,
type = true,
zero = false,
weight = 0
)
id = findall(x -> isa(x, Number), tbl.y)
tbl = tbl[id]
table = :id in propertynames(tbl) ? tbl : Table((; id = collect(1:length(tbl)), ), tbl)
table = :include in propertynames(tbl) ? table : Table(table; include = trues(length(table)))
f = get_formula(type, zero)
model = calfit(table, f, type, zero, weight)
xlevel = unique(table.x)
table = Table(; id = table.id, level = [findfirst(x -> i == x, xlevel) for i in table.x], y = table.y, x = table.x, x̂ = zeros(Float64, length(table)), accuracy = zeros(Float64, length(table)), include = table.include)
cal = Calibration(analyte, isd, type, zero, weight, f, source, table, model)
inv_predict_accuracy!(cal)
cal
end

function calfit(tbl, formula, type, zero, weight)
model = lm(formula, tbl[tbl.include]; wts = tbl.x[tbl.include] .^ weight)
if !type && !zero && model.model.pp.beta0[1] == 0
m = hcat(ones(Float64, count(tbl.include)), tbl.x[tbl.include], tbl.x[tbl.include] .^ 2)
sqrtw = diagm(sqrt.(tbl.x[tbl.include] .^ weight))
y = tbl.y[tbl.include]
model.model.pp.beta0 = (sqrtw * m) \ (sqrtw * y)
GLM.updateμ!(model.model.rr, predict(model, tbl[tbl.include]))
end
model
end

function calfit!(cal::MultipleCalibration)
cal.model = calfit(cal.table, cal.formula, cal.type, cal.zero, cal.weight)
cal
end

function project(cal::AbstractAnalysisTable, sample = nothing;
isd = 0,
type = true,
zero = false,
weight = 0
)
Project(
map(enumerate(cal.analytes)) do (id, analyte)
tbl = Table(;id = cal.samples, x = get_analyte(cal, analyte))
calibration(tbl, cal; analyte = id, isd, type, zero, weight)
end,
sample
)
end
72 changes: 72 additions & 0 deletions src/io.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@

function read_calibration(file::String, T, source::AbstractAnalysisTable, )
config_kw, config_vl = eachline(joinpath(file, "config.csv"))
config = Dict{Symbol, Any}()
for (k, v) in zip(split(config_kw, ","), split(config_vl, ","))
if k == "type" || k == "zero"
v = v == "TRUE" || v == "True" || v == "true" || v == ""
elseif k == "weight"
v = parse(Float64, v)
elseif k == "analyte" || k == "isd"
v = parse(Int, v)
elseif k == "beta"
v = parse(Float64, v)
else
continue
end
get!(config, Symbol(k), v)
end
endswith(file, ".scal") && return SingleCalibration(config[:analyte], config[:isd], source, config[:beta])
tbl = CSV.read(joinpath(file, "calibration.csv"), T)
calibration(tbl, source; config...)
end

function read_analysistable(file::String, T; anayte_fn = identity)
tbl = CSV.read(joinpath(file, "calibration.csv"), T)
config = filter!(!isempty, readlines(joinpath(file, "config.txt")))
config = replace.(config, Ref(" " => ""))
if pop!(config) == "R"
sample_name = Symbol(pop!(config))
analyte_name = Symbol.(config)
RowAnalysisTable(sample_name, analyte_name, analyte_fn.(analyte_name), tbl)
else
sample_name, analyte_name... = config
analyte_name = Symbol.(analyte_name)
ColumnAnalysisTable(Symbol(sample_name), analyte_name, analyte_fn.(analyte_name), tbl)
end
end

function read_project(file::String, T; analyte_fn = identity)
source = read_analysistable(file, T; analyte_fn)
fs = findfirst(f -> endswith(f, ".sample"), readdir(file))
Project([read_calibration(joinpath(file, f), T, source) for f in readdir(file) if endswith(f, ".cal") || endswith(f, ".scal")],
isnothing(fs) ? nothing : read_analysistable(joinpath(file, fs), T; analyte_fn)
)
end

function show(io::IO, cal::MultipleCalibration)
print(io, "Calibration of $(cal.source.analytes[cal.analyte]) with multiple levels")
end

function show(io::IO, cal::SingleCalibration)
print(io, "Calibration of $(cal.source.analytes[cal.analyte]) with single level")
end

function show(io::IO, ::MIME"text/plain", cal::MultipleCalibration),
print(io, cal, ":\n")
print(io, "∘ ISD: ", cal.isd > 0 ? cal.source.analytes[cal.isd] : nothing, "\n"),
print(io, "∘ Type: ", cal.type ? "quadratic" : "linear")
print(io, "∘ (0, 0): ", cal.zero ? "included\n" : "ommited\n")
print(io, "∘ Weight: ", weight_repr(cal.weight), "\n")
print(io, "∘ Formula", formula_repr(cal))
end

function show(io::IO, ::MIME"text/plain", cal::SingleCalibration)
print(io, cal, ":\n")
print(io, "∘ ISD: ", cal.isd > 0 ? cal.source.analytes[cal.isd] : nothing, "\n")
print(io, "∘ Formula", formula_repr(cal))
end

function show(io::io::IO, project::Project)
string("Project with $(length(project.calibration)) analytes")
end
91 changes: 91 additions & 0 deletions src/utils.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
get_formula(cal::MultipleCalibration) = get_formula(cal.type, cal.zero)
get_formula(type::Bool, zero::Bool) = if type
zero ? @formula(y ~ 0 + x) : @formula(y ~ x)
else
zero ? @formula(y ~ 0 + x + x ^ 2) : @formula(y ~ x + x ^ 2)
end

function get_analyte(tbl::RowAnalysisTable{T}, analyte::T) where T
id = findfirst(==(analyte), tbl.analytes)
isnothing(id) && throw(ArgumentError("Analyte $analyte is not in the table"))
collect(getproperties(tbl.table[id], tbl.sample_name))
end

function get_analyte(tbl::ColumnAnalysisTable{T}, analyte::T) where T
id = findfirst(==(analyte), tbl.analytes)
isnothing(id) && throw(ArgumentError("Analyte $analyte is not in the table"))
getproperty(tbl.table, tbl.analyte_name[id])
end

function get_sample(tbl::RowAnalysisTable{T}, sample::Symbol) where T
id = findfirst(==(sample), tbl.samples)
isnothing(id) && throw(ArgumentError("Sample $sample is not in the table"))
collect(getproperties(tbl.table[id], tbl.sample_name))
end

function get_sample(tbl::ColumnAnalysisTable{T}, sample::Symbol) where T
id = findfirst(==(sample), tbl.samples)
isnothing(id) && throw(ArgumentError("Sample $sample is not in the table"))
getproperty(tbl.table, tbl.analyte_name[id])
end

function critical_point(cal::MultipleCalibration)
β = cal.model.model.pp.beta0
c, b, a = cal.zero ? (0, β...) : β
-b / 2a
end

cal_range(project::Project) = cal_range(project.calibration)
cal_range(cal::MultipleCalibration) = (lloq(cal), uloq(cal))
lloq(project::Project) = lloq(project.calibration)
uloq(project::Project) = uloq(project.calibration)
lloq(cal::MultipleCalibration) = (cal.type || last(cal.model.model.pp.beta0) < 0) ? cal.table.x[findfirst(cal.table.include)] : max(cal.table.x[findfirst(cal.table.include)], critical_point(cal))
uloq(cal::MultipleCalibration) = (cal.type || last(cal.model.model.pp.beta0) > 0) ? cal.table.x[findlast(cal.table.include)] : min(cal.table.x[findlast(cal.table.include)], critical_point(cal))


function weight_repr(cal::MultipleCalibration)
cal.weight in [-0.5, -1, -2] || (cal.weight = 0)
weight_repr(cal.weight)
end
weight_repr(weight::Number) = if weight == -0.5
"1/√x"
elseif weight == -1
"1/x"
elseif weight == -2
"1/x²"
else
"none"
end

weight_value(weight) = if weight == "1/√x"
-0.5
elseif weight == "1/x"
-1
elseif weight == "1/x²"
-2
else
0
end

formula_repr(cal::SingleCalibration) = "y = $(round(; sigdigits = 4))x"
function formula_repr(cal::AbstractCalibration)
β = cal.model.model.pp.beta0
cal.type && cal.zero && return "y = $(round(β[1]; sigdigits = 4))x"
op = map(β[2:end]) do b
b < 0 ? " - " : " + "
end
if cal.type
string("y = ", format_number(β[1]), op[1], abs(format_number(β[2])), "x")
elseif cal.zero
string("y = ", format_number(β[1]), "x", op[1], abs(format_number(β[2])), "")
else
string("y = ", format_number(β[1]), op[1], abs(format_number(β[2])), "x", op[2], abs(format_number(β[3])), "")
end
end

formula_repr_utf8(cal::AbstractCalibration) = replace(formula_repr(cal), "" => "x^2")
weight_repr_utf8(cal::AbstractCalibration) = replace(weight_repr(cal), "" => "x^2", "√x" => "x^0.5")
format_number(x; digits) = format_number2int(round(x; digits))
format_number(x; sigdigits = 4) = format_number2int(round(x; sigdigits))
format_number2int(x) =
x == round(x) ? round(Int, x) : x

0 comments on commit bfa1d17

Please sign in to comment.