From ae8378cb23c4d697819671637be5145480eb7898 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Fri, 7 Feb 2020 15:35:55 +0300
Subject: [PATCH 01/15] plot update

---
 src/plots.jl | 111 ++++++++++++++++++++++++++++++---------------------
 1 file changed, 65 insertions(+), 46 deletions(-)

diff --git a/src/plots.jl b/src/plots.jl
index 573afdf..d3c807d 100644
--- a/src/plots.jl
+++ b/src/plots.jl
@@ -49,7 +49,7 @@ end
     xlabel      --> "Time"
     ylabel      --> "Concentration"
     link        --> :both
-    legend      --> false
+    legend      --> true
     grid        --> false
     #ticks       := [nothing :auto nothing]
     xlims       --> (minimum(x), maximum(x)),
@@ -68,12 +68,12 @@ function plotlabel(d)
         for (k, v) in d
             title *= "$(k) = $(v); "
         end
-        title = title[1:end-2]
+        title = title[1:length(title)]
     end
     return title
 end
 
-function plot(subj::T; title = nothing, legend = false, label = "AUTO", xlims = nothing, plotstyle::PKPlotStyle = PKPLOTSTYLE[1]) where T <: AbstractSubject
+function plot(subj::T; title = nothing, legend = true, label = "AUTO", xlims = nothing, plotstyle::PKPlotStyle = PKPLOTSTYLE[1]) where T <: AbstractSubject
 
     if title === nothing
         title = plotlabel(subj.sort)
@@ -92,7 +92,7 @@ function plot(subj::T; title = nothing, legend = false, label = "AUTO", xlims =
     return p
 
 end
-function plot!(subj::T; legend = false, label = "AUTO", xlims = nothing,  plotstyle::PKPlotStyle = PKPLOTSTYLE[1]) where T <: AbstractSubject
+function plot!(subj::T; legend = true, label = "AUTO", xlims = nothing,  plotstyle::PKPlotStyle = PKPLOTSTYLE[1]) where T <: AbstractSubject
     if xlims === nothing xlims = (minimum(subj.time), maximum(subj.time)) end
     p = pkplot!(subj.time, subj.obs;
         linestyle   = plotstyle.linestyle,
@@ -115,9 +115,54 @@ function usort(data::DataSet{T}, list) where T <: AbstractData
     dl
 end
 
-function plot(data::DataSet{T}; title = nothing, legend = false, xlims = nothing, pagesort = nothing, typesort = nothing) where T <: AbstractSubject
+function pageplot(pagedatatmp, styledict, typesort; title = "", legend = true, xlims = nothing)
+    utypes    = keys(styledict)
+    labels    = Vector{String}(undef, 0)
+    fst       = true
+    p         = nothing
+    for u in utypes
+        for d in pagedatatmp
+            if u ∈ d.sort
+                label = "AUTO"
+                if typesort !== nothing
+                    tempsd = Dict(k => d.sort[k] for k in typesort)
+                    style  = styledict[tempsd]
+                    label  = plotlabel(tempsd)
+                    if label ∉ labels
+                        push!(labels, label)
+                    else
+                        label = ""
+                    end
+                end
+                if fst
+                    p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = style, label = label)
+                    fst = false
+                else
+                    plot!(                  d; legend = legend, xlims = xlims, plotstyle = style, label = label)
+                end
+            end
+        end
+    end
+    p
+end
+function pageplot(pagedatatmp; title = "", legend = true, xlims = nothing)
+    fst       = true
+    p         = nothing
+    for d in pagedatatmp
+        label = "AUTO"
+        if fst
+            p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = PKPLOTSTYLE[1], label = label)
+            fst = false
+        else
+            plot!(                  d; legend = legend, xlims = xlims, plotstyle = PKPLOTSTYLE[1], label = label)
+        end
+    end
+    p
+end
+
+function plot(data::DataSet{T}; title = nothing, legend = true, xlims = nothing, pagesort = nothing, typesort = nothing) where T <: AbstractSubject
     # Style types
-    labels     = Vector{String}(undef, 0)
+
     styledict  = Dict()
     style      = PKPLOTSTYLE[1]
     if typesort !== nothing
@@ -133,59 +178,33 @@ function plot(data::DataSet{T}; title = nothing, legend = false, xlims = nothing
     autotitle = false
     if title === nothing autotitle = true end
     #---------------------------------------------------------------------------
+    # PAGES
     plots    = Vector{Any}(undef, 0)
     p        = nothing
     if pagesort !== nothing
-        pagedict = usort(data, pagesort)
+        pagedict   = usort(data, pagesort)
         for i in pagedict
             if autotitle
-                title = plotlabel(i)
+                title  = plotlabel(i)
             end
-            fst = true
+            pagedatatmp = Vector{eltype(data)}(undef, 0)
             for d in data
                 if i ∈ d.sort
-                    label = "AUTO"
-                    if typesort !== nothing
-                        tempsd = Dict(k => d.sort[k] for k in typesort)
-                        style  = styledict[tempsd]
-                        label = plotlabel(tempsd)
-                        if label ∉ labels
-                            push!(labels, label)
-                        else
-                            label = ""
-                        end
-                    end
-                    if fst
-                        p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = style, label = label)
-                        fst = false
-                    else
-                        plot!(d; plotstyle = style, legend = legend, xlims = xlims, label = label)
-                    end
+                    push!(pagedatatmp, d)
                 end
             end
-            push!(plots, p)
-        end
-    else
-        fst = true
-        for d in data
-            label = "AUTO"
             if typesort !== nothing
-                tempsd = Dict(k => d.sort[k] for k in typesort)
-                style  = styledict[tempsd]
-                label = plotlabel(tempsd)
-                if label ∉ labels
-                    push!(labels, label)
-                else
-                    label = ""
-                end
-            end
-
-            if fst
-                p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = style, label = label)
-                fst = false
+                p = pageplot(pagedatatmp, styledict, typesort; title = title, legend = legend, xlims = xlims)
             else
-                plot!(d; plotstyle = style, legend = legend, xlims = xlims, label = label)
+                p = pageplot(pagedatatmp; title = title, legend = legend, xlims = xlims)
             end
+            push!(plots, p)
+        end
+    else
+        if typesort !== nothing
+            p = pageplot(data, styledict, typesort; title = title, legend = legend, xlims = xlims)
+        else
+            p = pageplot(data; title = title, legend = legend, xlims = xlims)
         end
         push!(plots, p)
     end

From 5c457cfc1cddb9a45e71227437bd0889361c2ee0 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Fri, 7 Feb 2020 15:36:14 +0300
Subject: [PATCH 02/15] Cdose 0 if t0 0

---
 src/pk.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/pk.jl b/src/pk.jl
index db8a5d5..ae328a5 100644
--- a/src/pk.jl
+++ b/src/pk.jl
@@ -544,7 +544,12 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
                     result[:Cdose] = data.obs[i]
                     break
                 end
-                result[:Cdose] = linpredict(data.time[i] , data.time[i+1], data.dosetime.time, data.obs[i], data.obs[i+1])
+                #Cdose calculation
+                if data.obs[i] > 0
+                    result[:Cdose] = linpredict(data.time[i] , data.time[i+1], data.dosetime.time, data.obs[i], data.obs[i+1])
+                else
+                    result[:Cdose] = 0
+                end
                 aucpartl[i], aumcpartl[i] = aucpart(data.dosetime.time, data.time[i+1], result[:Cdose], data.obs[i+1], :lint, false) #? only :lint? always aftertmax = false?
                 break
             else

From 836288334fa09254774d79e99e30f8c68999dfd3 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Fri, 7 Feb 2020 15:36:50 +0300
Subject: [PATCH 03/15] eltype

---
 src/ClinicalTrialUtilities.jl | 1 +
 src/dataset.jl                | 3 +++
 2 files changed, 4 insertions(+)

diff --git a/src/ClinicalTrialUtilities.jl b/src/ClinicalTrialUtilities.jl
index d754e18..b43e10c 100644
--- a/src/ClinicalTrialUtilities.jl
+++ b/src/ClinicalTrialUtilities.jl
@@ -21,6 +21,7 @@ import Base.getindex
 import Base.length
 import Base.in
 import Base.iterate
+import Base.eltype
 import StatsBase.confint
 import DataFrames.DataFrame
 
diff --git a/src/dataset.jl b/src/dataset.jl
index 390c254..5a7ad7d 100644
--- a/src/dataset.jl
+++ b/src/dataset.jl
@@ -7,3 +7,6 @@ end
 function Base.iterate(iter::DataSet, i::Int)
     return Base.iterate(iter.data, i)
 end
+function Base.eltype(obj::DataSet)
+    return eltype(obj.data)
+end

From 71ed43edede60a3380f404c6b618c3699aa45ce4 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Sun, 9 Feb 2020 16:59:38 +0300
Subject: [PATCH 04/15] ConTab

---
 cange.log     |  5 ++++-
 src/freque.jl | 55 ++++++++++++++++++++++++++++++++++++++-------------
 test/test.jl  |  2 +-
 3 files changed, 46 insertions(+), 16 deletions(-)

diff --git a/cange.log b/cange.log
index da67ee6..8cd7bf9 100644
--- a/cange.log
+++ b/cange.log
@@ -1,3 +1,7 @@
+v0.2.3
+
+  - ConTab struct
+
 
 v0.2.2
   - Distributions bump
@@ -8,7 +12,6 @@ v0.2.2
   - fisher & pirson non-public methods for frequencies
   - add Show methods
   - minor fix
-  - cosmetics
   - add sim
   - cosmetics
 
diff --git a/src/freque.jl b/src/freque.jl
index 7ec50f4..d13e642 100644
--- a/src/freque.jl
+++ b/src/freque.jl
@@ -1,4 +1,26 @@
 
+
+struct ConTab{Int32, Int32}
+    tab::Matrix{Int}
+    row::Vector
+    col::Vector
+
+    function ConTab(m::Matrix{Int})
+        row = Vector{Symbol}(undef, size(m, 1)) .= Symbol("")
+        col = Vector{Symbol}(undef, size(m, 2)) .= Symbol("")
+        new{size(m, 1), size(m, 2)}(m, row, col)
+    end
+    function ConTab(m::Matrix{Int}, row, col)
+        new{size(m, 1), size(m, 2)}(m, row, col)
+    end
+end
+
+struct Freque
+    val
+    n
+end
+
+
 function freque(data::DataFrame; vars::Symbol, alpha = 0.05)::DataFrame
     result = DataFrame(value = Any[], n = Int[], p = Float64[], cil = Float64[], ciu = Float64[])
     list   = unique(data[:, vars])
@@ -12,12 +34,12 @@ function freque(data::DataFrame; vars::Symbol, alpha = 0.05)::DataFrame
     return result
 end
 
-function contab(data::DataFrame; row::Symbol, col::Symbol)::Matrix
+function contab(data::DataFrame; row::Symbol, col::Symbol)::ConTab
     clist = unique(data[:, col])
     rlist = unique(data[:, row])
     cn    = length(clist)
     rn    = length(rlist)
-    dfs   = Array{Int, 2}(undef, rn, cn)
+    dfs   = Matrix{Int}(undef, rn, cn)
     for ri = 1:rn
         rowl  = data[data[:, row] .== rlist[ri], col]
         for ci = 1:cn
@@ -25,29 +47,29 @@ function contab(data::DataFrame; row::Symbol, col::Symbol)::Matrix
             dfs[ri, ci] = cnt
         end
     end
-    return dfs
+    return ConTab(dfs, rlist, clist)
 end
 
 
 function pirson(a::Matrix{Int})
-    n = length(a[:,1])
-    m = length(a[1,:])
-    tm = sum(a, dims=1)[1,:]
-    tn = sum(a, dims=2)[:,1]
+    n   = length(a[:,1])
+    m   = length(a[1,:])
+    tm  = sum(a, dims=1)[1,:]
+    tn  = sum(a, dims=2)[:,1]
     num = sum(tm)
-    ae = Array{Real, 2}(undef, n, m)
+    ae  = Array{Real, 2}(undef, n, m)
     for im = 1:m
         for in = 1:n
             ae[in, im] = tn[in]*tm[im]/num
         end
     end
-    chsq = sum(((a .- ae) .^2 ) ./ ae)
+    chsq  = sum(((a .- ae) .^2 ) ./ ae)
     chsqy = sum(((abs.((a .- ae)) .- 0.5) .^2 ) ./ ae)
-    ml = 2 * sum( a .* log.( a ./ ae ))
-    df   = (n - 1)*(m - 1)
-    ϕ = sqrt(chsq / (num*(n - 1)*(m - 1)))
-    C = sqrt(chsq/(chsq+num))
-    K = sqrt(chsq/num/sqrt((n - 1)*(m - 1)))
+    ml    = 2 * sum( a .* log.( a ./ ae ))
+    df    = (n - 1)*(m - 1)
+    ϕ     = sqrt(chsq / (num*(n - 1)*(m - 1)))
+    C     = sqrt(chsq/(chsq+num))
+    K     = sqrt(chsq/num/sqrt((n - 1)*(m - 1)))
     return chsq, chsqy, ml, 1-cdf(Chisq(df), chsq)
 end
 
@@ -55,3 +77,8 @@ function fisher(a::Matrix{Int})
     dist  = Hypergeometric(sum(a[1, :]), sum(a[2, :]), sum(a[:, 1]))
     value = min(2 * min(cdf(dist, a[1, 1]), ccdf(dist, a[1, 1])), 1.0)
 end
+
+
+function fisher(t::ConTab{2, 2})
+    fisher(t.tab)
+end
diff --git a/test/test.jl b/test/test.jl
index b2de156..942e6b6 100644
--- a/test/test.jl
+++ b/test/test.jl
@@ -261,7 +261,7 @@ println(" ---------------------------------- ")
 @testset "  Frequency             " begin
     df   = freqdat
     ctab =  ClinicalTrialUtilities.contab(df, row = :row, col = :col)
-    @test ctab == [9 8; 5 21]
+    @test ctab.tab == [9 8; 5 21]
 
     frtab =  ClinicalTrialUtilities.freque(df; vars=:row, alpha = 0.05)
     @test frtab[1,2] == 17

From 98aa7395e2cb127f123a7e8944069ea4b9dc08a1 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Sun, 9 Feb 2020 17:01:05 +0300
Subject: [PATCH 05/15] types

---
 src/pk.jl | 50 ++++++++++++++++++++++++++++++++------------------
 1 file changed, 32 insertions(+), 18 deletions(-)

diff --git a/src/pk.jl b/src/pk.jl
index ae328a5..93f4fd6 100644
--- a/src/pk.jl
+++ b/src/pk.jl
@@ -76,9 +76,10 @@ struct DoseTime
 end
 
 
+#Plasma PK subject
 mutable struct PKSubject <: AbstractSubject
-    time::Array
-    obs::Array
+    time::Vector
+    obs::Vector
     kelauto::Bool
     kelrange::ElimRange
     dosetime::DoseTime
@@ -100,9 +101,22 @@ mutable struct PKSubject <: AbstractSubject
     end
 end
 
+#Urine PK subject
+mutable struct UPKSubject <: AbstractSubject
+    stime::Vector
+    etime::Vector
+    obs::Vector
+    vol::Vector
+    kelrange::ElimRange
+    dosetime::DoseTime
+    keldata::KelData
+    sort::Dict
+end
+
+#PD subject
 mutable struct PDSubject <: AbstractSubject
-    time::Array
-    obs::Array
+    time::Vector
+    obs::Vector
     bl::Real
     th::Real
     sort::Dict
@@ -369,7 +383,7 @@ end
     end
 
     #---------------------------------------------------------------------------
-    function slope(x::Vector, y::Vector)::Tuple{Number, Number, Number, Number}
+    function slope(x::Vector, y::Vector)
         if length(x) != length(y) throw(ArgumentError("Unequal vector length!")) end
         n   = length(x)
         if n < 2 throw(ArgumentError("n < 2!")) end
@@ -389,32 +403,32 @@ end
         return a, b, r2, ar
     end #end slope
         #Linear trapezoidal auc
-    function linauc(t₁, t₂, c₁, c₂)::Number
+    function linauc(t₁, t₂, c₁, c₂)
         return (t₂-t₁)*(c₁+c₂)/2
     end
         #Linear trapezoidal aumc
-    function linaumc(t₁, t₂, c₁, c₂)::Number
+    function linaumc(t₁, t₂, c₁, c₂)
         return (t₂-t₁)*(t₁*c₁+t₂*c₂)/2
     end
         #Log trapezoidal auc
-    function logauc(t₁, t₂, c₁, c₂)::Number
+    function logauc(t₁, t₂, c₁, c₂)
         return  (t₂-t₁)*(c₂-c₁)/log(c₂/c₁)
     end
         #Log trapezoidal aumc
-    function logaumc(t₁, t₂, c₁, c₂)::Number
+    function logaumc(t₁, t₂, c₁, c₂)
         return (t₂-t₁) * (t₂*c₂-t₁*c₁) / log(c₂/c₁) - (t₂-t₁)^2 * (c₂-c₁) / log(c₂/c₁)^2
     end
     #Intrapolation
         #linear prediction bx from ax, a1 < ax < a2
-    function linpredict(a₁, a₂, ax, b₁, b₂)::Number
+    function linpredict(a₁, a₂, ax, b₁, b₂)
         return abs((ax - a₁) / (a₂ - a₁))*(b₂ - b₁) + b₁
     end
 
-    function logtpredict(c₁, c₂, cx, t₁, t₂)::Number
+    function logtpredict(c₁, c₂, cx, t₁, t₂)
         return log(cx/c₁)/log(c₂/c₁)*(t₂-t₁)+t₁
     end
 
-    function logcpredict(t₁, t₂, tx, c₁, c₂)::Number
+    function logcpredict(t₁, t₂, tx, c₁, c₂)
         return exp(log(c₁) + abs((tx-t₁)/(t₂-t₁))*(log(c₂) - log(c₁)))
     end
     function cpredict(t₁, t₂, tx, c₁, c₂, calcm)
@@ -806,7 +820,7 @@ Pharmacokinetics data import from DataFrame.
 - conc - concentration column;
 - time - time column.
 """
-function pkimport(data::DataFrame, sort::Array, rule::LimitRule; conc::Symbol, time::Symbol)
+function pkimport(data::DataFrame, sort::Array, rule::LimitRule; conc::Symbol, time::Symbol)::DataSet
         sortlist = unique(data[:, sort])
         results  = Array{PKSubject, 1}(undef, 0)
         for i = 1:size(sortlist, 1) #For each line in sortlist
@@ -844,7 +858,7 @@ Pharmacokinetics data import from DataFrame.
 - conc - concentration column;
 - time - time column.
 """
-function pkimport(data::DataFrame, sort::Array; time::Symbol, conc::Symbol)
+function pkimport(data::DataFrame, sort::Array; time::Symbol, conc::Symbol)::DataSet
         rule = LimitRule()
         return pkimport(data, sort, rule; conc = conc, time = time)
 end
@@ -859,7 +873,7 @@ Pharmacokinetics data import from DataFrame for one subject.
 - conc - concentration column;
 - time - time column.
 """
-function pkimport(data::DataFrame, rule::LimitRule; time::Symbol, conc::Symbol)
+function pkimport(data::DataFrame, rule::LimitRule; time::Symbol, conc::Symbol)::DataSet
         #rule = LimitRule()
         datai = ncarule!(copy(data[!,[time, conc]]), conc, time, rule)
         return DataSet([PKSubject(datai[!, time], datai[!, conc])])
@@ -874,12 +888,12 @@ Pharmacokinetics data import from DataFrame for one subject.
 - conc - concentration column;
 - time - time column.
 """
-function pkimport(data::DataFrame; time::Symbol, conc::Symbol)
+function pkimport(data::DataFrame; time::Symbol, conc::Symbol)::DataSet
         datai = sort(data[!,[time, conc]], time)
         return DataSet([PKSubject(datai[!, time], datai[!, conc])])
 end
 
-function pkimport(time::Vector, conc::Vector; sort = Dict())
+function pkimport(time::Vector, conc::Vector; sort = Dict())::PKSubject
      PKSubject(time, conc; sort = sort)
 end
     #---------------------------------------------------------------------------
@@ -897,7 +911,7 @@ Pharmacodynamics data import from DataFrame.
 - bl - baseline;
 - th - treashold.
 """
-function pdimport(data::DataFrame, sort::Array; resp::Symbol, time::Symbol, bl::Real = 0, th::Real = NaN)
+function pdimport(data::DataFrame, sort::Array; resp::Symbol, time::Symbol, bl::Real = 0, th::Real = NaN)::DataSet
         sortlist = unique(data[:, sort])
         results  = Array{PDSubject, 1}(undef, 0)
         for i = 1:size(sortlist, 1) #For each line in sortlist

From 33428da91f47e1473ff76398b008c5e09485edd6 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Tue, 11 Feb 2020 21:38:35 +0300
Subject: [PATCH 06/15] DataFrames 0.20

---
 Project.toml | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 06617d8..6a18237 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,6 +12,7 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [extras]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -25,9 +26,10 @@ test = ["CSV", "Test", "Plots"]
 julia = "1.0, 1.1, 1.2, 1.3"
 Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22"
 StatsBase = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32"
-DataFrames = "0.19"
+DataFrames = "0.19, 0.20"
 QuadGK = "2.0, 2.1"
 SpecialFunctions = "0.8, 0.9"
 Roots = "0.7, 0.8"
 CSV = "0.5"
 RecipesBase = "0.7"
+Reexport = "0.1, 0.2"

From b386a27454faffdbda757b0e1eba0088f60245c5 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Tue, 11 Feb 2020 21:39:17 +0300
Subject: [PATCH 07/15] minor fix

---
 src/ClinicalTrialUtilities.jl |  9 ++++++---
 src/descriptives.jl           | 10 +++++++++-
 src/plots.jl                  | 11 ++++++++---
 test/pktest.jl                |  4 ++--
 test/test.jl                  |  2 +-
 5 files changed, 26 insertions(+), 10 deletions(-)

diff --git a/src/ClinicalTrialUtilities.jl b/src/ClinicalTrialUtilities.jl
index b43e10c..3a43ad3 100644
--- a/src/ClinicalTrialUtilities.jl
+++ b/src/ClinicalTrialUtilities.jl
@@ -13,7 +13,11 @@
 # If you want to check and get R code you can find some here: http://powerandsamplesize.com/Calculators/
 __precompile__(true)
 module ClinicalTrialUtilities
-using Distributions, StatsBase, Random, Roots, QuadGK, DataFrames, RecipesBase
+
+using Distributions, Random, Roots, QuadGK, RecipesBase, Reexport
+
+@reexport using StatsBase
+
 import SpecialFunctions
 import Base.show
 import Base.showerror
@@ -23,8 +27,7 @@ import Base.in
 import Base.iterate
 import Base.eltype
 import StatsBase.confint
-import DataFrames.DataFrame
-
+import DataFrames: DataFrame, DataFrames, names!, unstack, deleterows!
 
 function lgamma(x)
     return SpecialFunctions.logabsgamma(x)[1]
diff --git a/src/descriptives.jl b/src/descriptives.jl
index 61cef2b..6953752 100644
--- a/src/descriptives.jl
+++ b/src/descriptives.jl
@@ -183,13 +183,21 @@ end
 
 #Statistics calculation
 
-@inline function descriptive_(data::Array{T,1}, stats::Union{Tuple{Vararg{Symbol}}, Array{Symbol,1}}) where T <: Real
+function notnan(x)
+    return !(x === NaN || x === nothing || x === missing)
+end
+
+@inline function descriptive_(data::Vector{T}, stats::Union{Tuple{Vararg{Symbol}}, Array{Symbol,1}}) where T <: Real
 
+    #=
     dlist = findall(x -> x === NaN || x === nothing || x === missing, data)
     if length(dlist) > 0
         data = copy(data)
         deleteat!(data, dlist)
     end
+    =#
+
+    data       = data[notnan.(data)]
 
     dn         = nothing
     dmin       = nothing
diff --git a/src/plots.jl b/src/plots.jl
index d3c807d..6db1ffa 100644
--- a/src/plots.jl
+++ b/src/plots.jl
@@ -56,7 +56,7 @@ end
     ylims       --> (0, maximum(y)*1.1)
     seriescolor --> :blue
     markershape --> :circle
-    markersize  --> 2
+    markersize  --> 3
     markercolor --> :match
     msalpha     --> 0
     (x, y)
@@ -163,9 +163,14 @@ end
 function plot(data::DataSet{T}; title = nothing, legend = true, xlims = nothing, pagesort = nothing, typesort = nothing) where T <: AbstractSubject
     # Style types
 
-    styledict  = Dict()
-    style      = PKPLOTSTYLE[1]
+    styledict  = nothing
+    #style      = PKPLOTSTYLE[1]
+    if pagesort !== nothing
+        if !isa(pagesort, Array) pagesort = [pagesort] end
+    end
     if typesort !== nothing
+        styledict  = Dict()
+        if !isa(typesort, Array) typesort = [typesort] end
         utypes = usort(data, typesort)
         for i = 1:length(utypes)
             if i <= 20
diff --git a/test/pktest.jl b/test/pktest.jl
index 0c31ee5..072d75d 100644
--- a/test/pktest.jl
+++ b/test/pktest.jl
@@ -356,8 +356,8 @@ println(" ---------------------------------- ")
     pkds = ClinicalTrialUtilities.pkimport(glucose2, [:Subject, :Date]; conc = :glucose, time = :Time)
     pk   = ClinicalTrialUtilities.nca!(pkds)
     df   = DataFrame(pk; unst = true)
-
-
+    p    = ClinicalTrialUtilities.plot(pkds; pagesort = [:Date], typesort = [:Subject])
+    @test length(p) == 2
     print(io, pk[1].subject.keldata)
 end
 
diff --git a/test/test.jl b/test/test.jl
index 942e6b6..68bc360 100644
--- a/test/test.jl
+++ b/test/test.jl
@@ -1,6 +1,6 @@
 # Clinical Trial Utilities
 # Copyright © 2019 Vladimir Arnautov aka PharmCat (mail@pharmcat.net)
-using Distributions, Random, DataFrames, CSV, Test
+using Distributions, Random, DataFrames, CSV, Test, Plots
 
 path    = dirname(@__FILE__)
 io      = IOBuffer();

From 3d79fc9c2a17027579ded23707737bc987894548 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Tue, 11 Feb 2020 21:40:28 +0300
Subject: [PATCH 08/15] minor improve, verbose

---
 src/pk.jl | 277 +++++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 203 insertions(+), 74 deletions(-)

diff --git a/src/pk.jl b/src/pk.jl
index 93f4fd6..ba02c0d 100644
--- a/src/pk.jl
+++ b/src/pk.jl
@@ -522,6 +522,9 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
         result[:Obsnum]  = length(data)
         result[:Cmax]    = maximum(data.obs)
 
+        doseaucpart  = 0.0
+        doseaumcpart = 0.0
+
         for i = result[:Obsnum]:-1:1
             if data.obs[i] > 0
                 result[:Tlast]   = data.time[i]
@@ -560,34 +563,24 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
                 end
                 #Cdose calculation
                 if data.obs[i] > 0
-                    result[:Cdose] = linpredict(data.time[i] , data.time[i+1], data.dosetime.time, data.obs[i], data.obs[i+1])
+                    result[:Cdose] = linpredict(data.time[i] , data.time[i + 1], data.dosetime.time, data.obs[i], data.obs[i + 1])
                 else
                     result[:Cdose] = 0
                 end
-                aucpartl[i], aumcpartl[i] = aucpart(data.dosetime.time, data.time[i+1], result[:Cdose], data.obs[i+1], :lint, false) #? only :lint? always aftertmax = false?
+                #aucpartl[i], aumcpartl[i] = aucpart(data.dosetime.time, data.time[i+1], result[:Cdose], data.obs[i+1], :lint, false) #? only :lint? always aftertmax = false?
+                doseaucpart, doseaumcpart  = aucpart(data.dosetime.time, data.time[i + 1], result[:Cdose], data.obs[i + 1], :lint, false)
+                pmask[i]     = false
                 break
             else
-                aucpartl[i]  = 0
-                aumcpartl[i] = 0
+                #aucpartl[i]  = 0
+                #aumcpartl[i] = 0
                 pmask[i]     = false
             end
         end
         #sum all AUC parts
-        result[:AUCall]  = sum(aucpartl[pmask])
-        result[:AUMCall] = sum(aumcpartl[pmask])
+        result[:AUCall]  = sum(aucpartl[pmask])  + doseaucpart
+        result[:AUMCall] = sum(aumcpartl[pmask]) + doseaumcpart
         #-----------------------------------------------------------------------
-        if verbose
-            aucpartlsum  = similar(aucpartl)
-            aumcpartlsum = similar(aumcpartl)
-            for i = 1:length(aucpartl)
-                aucpartlsum[i]  = sum(aucpartl[1:i])
-                aumcpartlsum[i] = sum(aumcpartl[1:i])
-            end
-            mx = hcat(data.time, data.obs, round.(vcat([0], aucpartl), digits = 3),  round.(vcat([0], aucpartlsum), digits = 3), round.(vcat([0], aumcpartl), digits = 3),  round.(vcat([0], aumcpartlsum), digits = 3))
-            mx = vcat(permutedims(["Time", "Concentrtion", "AUC", "AUC (cumulate)", "AUMC", "AUMC (cumulate)"]), mx)
-            printmatrix(io, mx)
-            println(io, "")
-        end
         #-----------------------------------------------------------------------
         #Exclude all AUC parts from observed concentation before 0 or less
         #Need elaborate!!!!
@@ -641,23 +634,19 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
         end
         #-----------------------------------------------------------------------
         #Steady-state
+        tautime    = data.dosetime.time + data.dosetime.tau
         if data.dosetime.tau > 0
             ncas, ncae = ncarange(data, data.dosetime.time, data.dosetime.tau)
-            tautime    = data.dosetime.time + data.dosetime.tau
+            #tautime    = data.dosetime.time + data.dosetime.tau
             eaucpartl  = eaumcpartl = 0.0
             if tautime < data.time[end]
                 #result[:Ctau] = linpredict(data.time[ncae] , data.time[ncae+1], tautime, data.obs[ncae], data.obs[ncae+1])
                 if tautime > result[:Tmax] aftertmax = true else aftertmax = false end
-                result[:Ctau] = cpredict(data.time[ncae], data.time[ncae+1], tautime, data.obs[ncae], data.obs[ncae+1], intp)
-                aucpartl[ncae], aumcpartl[ncae] = aucpart(data.time[ncae], tautime, data.obs[ncae], result[:Ctau], calcm, aftertmax)
-                if verbose
-                    println(io, "Tau + dosetime is less then end time. Inrapolation used.")
-                    println(io, "Interpolation between: $(data.time[ncae]) - $( data.time[ncae+1]), method: $(intp)")
-                    println(io, "Ctau: $(result[:Ctau])")
-                    println(io, "AUC final part: $(aucpartl[ncae])")
-                end
+                result[:Ctau] = cpredict(data.time[ncae], data.time[ncae + 1], tautime, data.obs[ncae], data.obs[ncae + 1], intp)
+                eaucpartl, eaumcpartl = aucpart(data.time[ncae], tautime, data.obs[ncae], result[:Ctau], calcm, aftertmax)
+
                 #remoove part after tau
-                if ncae < result[:Obsnum] - 1 pmask[ncae+1:end] .= false end
+                if ncae < result[:Obsnum] - 1 pmask[ncae:end] .= false end
             elseif tautime > data.time[end]
                 #extrapolation
                 result[:Ctau] = exp(result[:LZint] + result[:LZ] * tautime)
@@ -685,6 +674,36 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
             result[:Fluctau]  = (result[:Cmax] - result[:Ctau])/result[:Cavg]
             #result[:AccInd]
         end
+        if verbose
+            aucpartlsum  = similar(aucpartl)
+            aumcpartlsum = similar(aumcpartl)
+            for i = 1:length(aucpartl)
+                aucpartlsum[i]  = sum(aucpartl[1:i])
+                aumcpartlsum[i] = sum(aumcpartl[1:i])
+            end
+            astx    = Vector{String}(undef, length(data.time))
+            astx[1] = ""
+            for i = 1:length(pmask)
+                if pmask[i] astx[i+1] = "Yes" else astx[i+1] = "No" end
+            end
+            mx = hcat(data.time, data.obs, round.(vcat([0], aucpartl), digits = 3),  round.(vcat([0], aucpartlsum), digits = 3), round.(vcat([0], aumcpartl), digits = 3),  round.(vcat([0], aumcpartlsum), digits = 3), astx)
+            mx = vcat(permutedims(["Time", "Concentrtion", "AUC", "AUC (cumulate)", "AUMC", "AUMC (cumulate)", "Include"]), mx)
+            printmatrix(io, mx)
+            println(io, "")
+            println(io, "Cdose: $(result[:Cdose]), Dose time: $(data.dosetime.time)")
+            if data.dosetime.time > data.time[1]
+                println("Dose AUC  part: $(doseaucpart)")
+                println("Dose AUMC part: $(doseaumcpart)")
+            end
+            println(io, "")
+            if tautime < data.time[end]
+                println(io, "Tau + dosetime is less then end time. Interpolation used.")
+                println(io, "Interpolation between: $(data.time[ncae]) - $( data.time[ncae+1]), method: $(intp)")
+                println(io, "Ctau: $(result[:Ctau])")
+                println(io, "AUC  final part: $(eaucpartl)")
+                println(io, "AUMC final part: $(eaumcpartl)")
+            end
+        end
         #-----------------------------------------------------------------------
         return PKPDProfile(data, result; method = calcm)
     end
@@ -806,8 +825,59 @@ function nca!(data::DataSet{PDSubject}; verbose = false, io::IO = stdout)
             push!(results, nca!(data.data[i]))
         end
         return DataSet(results)
+end
+
+
+"""
+    nca!(data::PDSubject; verbose = false, io::IO = stdout)::PKPDProfile{PDSubject}
+
+Pharmacodynamics non-compartment analysis for urine data.
+"""
+function nca!(data::UPKSubject; verbose = false, io::IO = stdout)::PKPDProfile{UPKSubject}
+
+    result   = Dict()
+    result[:midtime] = (data.stime .+ data.etime) / 2
+    result[:exrate]  = (data.obs .* data.vol) ./  (data.stime .+ data.etime)
+    result[:ar]      =  sum(data.obs .* data.vol)
+    result[:maxrate] = maximum(result[:exrate])
+    result[:volume]  = sum(data.vol)
+    #result[:AUCrate]
+    #result[:tmaxrate]
+    #result[:arp]
+    #result[:Kel]
+    #result[:HL]
+end
+#=
+stime::Vector
+etime::Vector
+obs::Vector
+vol::Vector
+=#
+#-------------------------------------------------------------------------------
+
+
+function getdatai(data, sort, cols, func; sortby = nothing)
+    sortlist = unique(data[:, sort])
+    datai    = DataFrame()
+    for i in cols
+        datai[!, i] = Vector{eltype(data[!, i])}(undef, 0)
     end
-    #---------------------------------------------------------------------------
+    for i = 1:size(sortlist, 1)
+        if size(datai, 1) > 0 deleterows!(datai, 1:size(datai, 1)) end
+        for c = 1:size(data, 1) #For each line in data
+            if data[c, sort] == sortlist[i,:]
+                push!(datai, data[c, cols])
+            end
+        end
+        if sortby !== nothing
+            sort!(datai, sortby)
+        end
+        func(datai, sortlist[i,:])
+    end
+end
+
+#-------------------------------------------------------------------------------
+
 """
     pkimport(data::DataFrame, sort::Array, rule::LimitRule; conc::Symbol, time::Symbol)
 
@@ -821,32 +891,39 @@ Pharmacokinetics data import from DataFrame.
 - time - time column.
 """
 function pkimport(data::DataFrame, sort::Array, rule::LimitRule; conc::Symbol, time::Symbol)::DataSet
-        sortlist = unique(data[:, sort])
-        results  = Array{PKSubject, 1}(undef, 0)
-        for i = 1:size(sortlist, 1) #For each line in sortlist
-            datai = DataFrame(Time = eltype(data[!,time])[], Conc = eltype(data[!,conc])[])
-            for c = 1:size(data, 1) #For each line in data
-                if data[c, sort] == sortlist[i,:]
-                    push!(datai, [data[c, time], data[c, conc]])
-                end
-            end
-            if !(eltype(datai.Time) <: Real)
-                datai.Time[.!isa.(datai.Time, Real)] .= NaN
-                datai.Time = float.(datai.Time)
-            end
-            if !(eltype(datai.Conc) <: Real)
-                datai.Conc[.!isa.(datai.Conc, Real)] .= NaN
-                datai.Conc = float.(datai.Conc)
+    results  = Array{PKSubject, 1}(undef, 0)
+    getdatai(data, sort, [conc, time], (x, y) -> push!(results, PKSubject(copy(x[!, time]), copy(x[!, conc]), sort = Dict(sort .=> collect(y)))); sortby = time)
+    #=
+    sortlist = unique(data[:, sort])
+    for i = 1:size(sortlist, 1) #For each line in sortlist
+        datai = DataFrame(Time = eltype(data[!,time])[], Conc = eltype(data[!,conc])[])
+        for c = 1:size(data, 1) #For each line in data
+            if data[c, sort] == sortlist[i,:]
+                push!(datai, [data[c, time], data[c, conc]])
             end
-            if rule.lloq !== NaN
-                ncarule!(datai, :Conc, :Time, rule)
-            else
-                sort!(datai, :Time)
-            end
-            push!(results, PKSubject(datai[!, :Time], datai[!, :Conc], sort = Dict(sort .=> collect(sortlist[i,:]))))
         end
-        return DataSet(results)
+        if !(eltype(datai.Time) <: Real)
+            datai.Time[.!isa.(datai.Time, Real)] .= NaN
+            datai.Time = float.(datai.Time)
+        end
+        if !(eltype(datai.Conc) <: Real)
+            datai.Conc[.!isa.(datai.Conc, Real)] .= NaN
+            datai.Conc = float.(datai.Conc)
+        end
+        if rule.lloq !== NaN
+            ncarule!(datai, :Conc, :Time, rule)
+        else
+            sort!(datai, :Time)
+        end
+        push!(results, PKSubject(datai[!, :Time], datai[!, :Conc], sort = Dict(sort .=> collect(sortlist[i,:]))))
+    end
+    =#
+    results = DataSet(results)
+    if rule.lloq !== NaN
+         applyncarule!(results, rule)
     end
+    return results
+end
 """
     pkimport(data::DataFrame, sort::Array; time::Symbol, conc::Symbol)
 
@@ -859,8 +936,8 @@ Pharmacokinetics data import from DataFrame.
 - time - time column.
 """
 function pkimport(data::DataFrame, sort::Array; time::Symbol, conc::Symbol)::DataSet
-        rule = LimitRule()
-        return pkimport(data, sort, rule; conc = conc, time = time)
+    rule = LimitRule()
+    return pkimport(data, sort, rule; conc = conc, time = time)
 end
 """
     pkimport(data::DataFrame, sort::Array; time::Symbol, conc::Symbol)
@@ -875,8 +952,8 @@ Pharmacokinetics data import from DataFrame for one subject.
 """
 function pkimport(data::DataFrame, rule::LimitRule; time::Symbol, conc::Symbol)::DataSet
         #rule = LimitRule()
-        datai = ncarule!(copy(data[!,[time, conc]]), conc, time, rule)
-        return DataSet([PKSubject(datai[!, time], datai[!, conc])])
+    datai = ncarule!(copy(data[!,[time, conc]]), conc, time, rule)
+    return DataSet([PKSubject(datai[!, time], datai[!, conc])])
 end
 """
     pkimport(data::DataFrame, sort::Array; time::Symbol, conc::Symbol)
@@ -889,8 +966,8 @@ Pharmacokinetics data import from DataFrame for one subject.
 - time - time column.
 """
 function pkimport(data::DataFrame; time::Symbol, conc::Symbol)::DataSet
-        datai = sort(data[!,[time, conc]], time)
-        return DataSet([PKSubject(datai[!, time], datai[!, conc])])
+    datai = sort(data[!,[time, conc]], time)
+    return DataSet([PKSubject(datai[!, time], datai[!, conc])])
 end
 
 function pkimport(time::Vector, conc::Vector; sort = Dict())::PKSubject
@@ -912,41 +989,82 @@ Pharmacodynamics data import from DataFrame.
 - th - treashold.
 """
 function pdimport(data::DataFrame, sort::Array; resp::Symbol, time::Symbol, bl::Real = 0, th::Real = NaN)::DataSet
-        sortlist = unique(data[:, sort])
-        results  = Array{PDSubject, 1}(undef, 0)
-        for i = 1:size(sortlist, 1) #For each line in sortlist
-            datai = DataFrame(Time = Real[], Resp = Real[])
-            for c = 1:size(data, 1) #For each line in data
-                if data[c, sort] == sortlist[i,:]
-                    push!(datai, [data[c, time], data[c, resp]])
-                end
+    №sortlist = unique(data[:, sort])
+    results  = Array{PDSubject, 1}(undef, 0)
+    getdatai(data, sort, [resp, time], (x, y) -> push!(results, PDSubject(copy(x[!, time]), copy(x[!, resp]), bl, th, sort = Dict(sort .=> collect(y)))); sortby = time)
+    #=
+    for i = 1:size(sortlist, 1) #For each line in sortlist
+        datai = DataFrame(Time = Real[], Resp = Real[])
+        for c = 1:size(data, 1) #For each line in data
+            if data[c, sort] == sortlist[i,:]
+                push!(datai, [data[c, time], data[c, resp]])
             end
-            sort!(datai, :Time)
-            push!(results, PDSubject(datai[!, :Time], datai[!, :Resp], bl, th, sort = Dict(sort .=> collect(sortlist[i,:]))))
         end
-        return DataSet(results)
-end
-    function pdimport(data::DataFrame; resp::Symbol, time::Symbol, bl = 0, th = NaN)
-        datai = sort(data[!,[time, resp]], time)
-        return DataSet([PDSubject(datai[!, time], datai[!, resp], bl, th)])
+        sort!(datai, :Time)
+        push!(results, PDSubject(datai[!, :Time], datai[!, :Resp], bl, th, sort = Dict(sort .=> collect(sortlist[i,:]))))
     end
-    #---------------------------------------------------------------------------
+    =#
+    return DataSet(results)
+end
+
+function pdimport(data::DataFrame; resp::Symbol, time::Symbol, bl = 0, th = NaN)
+    datai = sort(data[!,[time, resp]], time)
+    return DataSet([PDSubject(datai[!, time], datai[!, resp], bl, th)])
+end
+
+#-------------------------------------------------------------------------------
     function applyncarule(data::PKSubject, rule::LimitRule)
+        cmax, tmax, tmaxn = ctmax(data)
+        #NaN Rule
+        obsn = length(data)
+        if rule.nan !== NaN
+            for i = 1:length(data)
+                if data.obs[i] === NaN data.obs[i] = rule.nan end
+            end
+        end
+
+        #LLOQ rule
+        for i = 1:obsn
+            if data.obs[i] < rule.lloq
+                if i <= tmaxn
+                    data.obs[i] = rule.btmax
+                else
+                    data.obs[i] = rule.atmax
+                end
+            end
+        end
+
+        #NaN Remove rule
+        if rule.rm
+            filterv   = data.obs .!== NaN
+            data.time = data.time[filterv]
+            data.obs  = data.obs[filterv]
+        end
     end
     function applyncarule(data::DataSet{PKSubject}, rule::LimitRule)
     end
+    #=
     function applyncarule(data::PKPDProfile{PKSubject}, rule::LimitRule)
     end
     function applyncarule(data::DataSet{PKPDProfile{PKSubject}}, rule::LimitRule)
     end
+    =#
+
     function applyncarule!(data::PKSubject, rule::LimitRule)
+
     end
     function applyncarule!(data::DataSet{PKSubject}, rule::LimitRule)
+        for i in data
+            applyncarule!(i, rule)
+        end
+        data
     end
+    #=
     function applyncarule!(data::PKPDProfile{PKSubject}, rule::LimitRule)
     end
     function applyncarule!(data::DataSet{PKPDProfile{PKSubject}}, rule::LimitRule)
     end
+    =#
     #---------------------------------------------------------------------------
 """
     setelimrange!(data::PKSubject, range::ElimRange)
@@ -1168,3 +1286,14 @@ function DataFrames.DataFrame(data::DataSet{T}; unst = false, us = false) where
             return df
         end
 end
+
+
+#-------------------------------------------------------------------------------
+# Param functions
+
+function cmax(data::PKSubject)
+    return maximum(data.obs)
+end
+
+function auc()
+end

From 34d664057ad629c758caffbba2cf0ec77cba531f Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Wed, 12 Feb 2020 22:03:30 +0300
Subject: [PATCH 09/15] PK fix and update

---
 src/pk.jl            | 167 +++++++++++++++++++++++++++----------------
 src/plots.jl         |  22 +++---
 test/csv/nanpk.csv   |  22 ++++++
 test/csv/upkdata.csv |   6 ++
 test/pktest.jl       |  45 +++++++++---
 test/testdata.jl     |   4 ++
 6 files changed, 186 insertions(+), 80 deletions(-)
 create mode 100644 test/csv/nanpk.csv
 create mode 100644 test/csv/upkdata.csv

diff --git a/src/pk.jl b/src/pk.jl
index ba02c0d..e87db60 100644
--- a/src/pk.jl
+++ b/src/pk.jl
@@ -105,12 +105,18 @@ end
 mutable struct UPKSubject <: AbstractSubject
     stime::Vector
     etime::Vector
-    obs::Vector
+    conc::Vector
     vol::Vector
+    time::Vector
+    obs::Vector
+    kelauto::Bool
     kelrange::ElimRange
     dosetime::DoseTime
     keldata::KelData
     sort::Dict
+    function UPKSubject(stime::Vector, etime::Vector, conc::Vector, vol::Vector; sort = Dict())
+        new(stime, etime, conc, vol, (stime .+ etime) / 2, (conc .* vol) ./  (etime .- stime), true, ElimRange(), DoseTime(NaN, 0 * stime[1]), KelData(), sort)::UPKSubject
+    end
 end
 
 #PD subject
@@ -159,7 +165,7 @@ struct LimitRule
 end
 
 
-function Base.getindex(a::PKPDProfile{T}, s::Symbol)::Real where T <: AbstractSubject
+function Base.getindex(a::PKPDProfile{T}, s::Symbol) where T <: AbstractSubject
     return a.result[s]
 end
 
@@ -300,6 +306,9 @@ function ncarule!(data::PKSubject, rule::LimitRule)
     #!!!!
 end
 
+function notnanormissing(x)
+    if x === NaN || x === missing return false else true end
+end
     """
         cmax
         tmax
@@ -307,14 +316,16 @@ end
     """
     function ctmax(time::Vector, obs::Vector)
         if length(time) != length(obs) throw(ArgumentError("Unequal vector length!")) end
-        if length(obs[obs .!== NaN]) > 0 cmax  = maximum(obs[obs .!== NaN]) else cmax = NaN end
+        if length(obs[notnanormissing.(obs)]) > 0 cmax  = maximum(obs[notnanormissing.(obs)]) else cmax = NaN end
         tmax  = NaN
         tmaxn = 0
         for i = 1:length(time)
-            if obs[i] == cmax
-                tmax  = time[i]
-                tmaxn = i
-                break
+            if obs[i] !== missing
+                if obs[i] == cmax
+                    tmax  = time[i]
+                    tmaxn = i
+                    break
+                end
             end
         end
         return cmax, tmax, tmaxn
@@ -367,10 +378,12 @@ end
     """
     function ncarange(data::PKSubject, dosetime, tau)
         tautime = dosetime + tau
-        s     = 0
-        e     = 0
+        s     = 1
+        e     = length(data.time)
         for i = 1:length(data) - 1
-            if dosetime >= data.time[i] && dosetime < data.time[i+1] s = i; break end
+            if dosetime > data.time[i] && dosetime <= data.time[i+1] s = i+1;
+                break
+            end
         end
         if tautime >= data.time[end]
             e = length(data.time)
@@ -551,46 +564,47 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
         pmask   = Array{Bool, 1}(undef, result[:Obsnum] - 1)
         pmask  .= true
 
+        #sum all AUC parts
+        result[:AUCall]  = sum(aucpartl)
+        result[:AUMCall] = sum(aumcpartl)
+        #-----------------------------------------------------------------------
+        #-----------------------------------------------------------------------
+        #Exclude all AUC parts from observed concentation before 0 or less
+        #Need elaborate!!!!
+        for i = result[:Tmaxn]:result[:Obsnum] - 1
+            if data.obs[i+1] <= 0 * data.obs[i+1] pmask[i:end] .= false; break end
+        end
+        result[:AUClast]  = sum(aucpartl[pmask])
+        result[:AUMClast] = sum(aumcpartl[pmask])
+        result[:MRTlast]       = result[:AUMClast] / result[:AUClast]
+
+        #-----------------------------------------------------------------------
+
         if data.dosetime.time < data.time[1]
             #Dosetime < first time
         end
         #Find AUC part from dose time; exclude all before dosetime
         for i = 1:result[:Obsnum] - 1
-            if  data.time[i] <= data.dosetime.time < data.time[i+1]
+            if  data.time[i] <= data.dosetime.time <= data.time[i+1]
                 if data.time[i] == data.dosetime.time
                     result[:Cdose] = data.obs[i]
                     break
-                end
-                #Cdose calculation
-                if data.obs[i] > 0
-                    result[:Cdose] = linpredict(data.time[i] , data.time[i + 1], data.dosetime.time, data.obs[i], data.obs[i + 1])
                 else
-                    result[:Cdose] = 0
-                end
+                #Cdose calculation
+                    if data.obs[i] > 0
+                        result[:Cdose] = linpredict(data.time[i] , data.time[i + 1], data.dosetime.time, data.obs[i], data.obs[i + 1])
+                    else
+                        result[:Cdose] = 0
+                    end
                 #aucpartl[i], aumcpartl[i] = aucpart(data.dosetime.time, data.time[i+1], result[:Cdose], data.obs[i+1], :lint, false) #? only :lint? always aftertmax = false?
-                doseaucpart, doseaumcpart  = aucpart(data.dosetime.time, data.time[i + 1], result[:Cdose], data.obs[i + 1], :lint, false)
-                pmask[i]     = false
-                break
+                    doseaucpart, doseaumcpart  = aucpart(data.dosetime.time, data.time[i + 1], result[:Cdose], data.obs[i + 1], :lint, false)
+                    pmask[i]     = false
+                    break
+                end
             else
-                #aucpartl[i]  = 0
-                #aumcpartl[i] = 0
                 pmask[i]     = false
             end
         end
-        #sum all AUC parts
-        result[:AUCall]  = sum(aucpartl[pmask])  + doseaucpart
-        result[:AUMCall] = sum(aumcpartl[pmask]) + doseaumcpart
-        #-----------------------------------------------------------------------
-        #-----------------------------------------------------------------------
-        #Exclude all AUC parts from observed concentation before 0 or less
-        #Need elaborate!!!!
-        for i = result[:Tmaxn]:result[:Obsnum] - 1
-            if data.obs[i+1] <= 0 * data.obs[i+1] pmask[i:end] .= false; break end
-        end
-        result[:AUClast]  = sum(aucpartl[pmask])
-        result[:AUMClast] = sum(aumcpartl[pmask])
-
-        result[:MRTlast]       = result[:AUMClast] / result[:AUClast]
         #-----------------------------------------------------------------------
         #Elimination
         keldata                = KelData()
@@ -630,6 +644,7 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
             result[:Clast_pred]      = exp(result[:LZint] + result[:LZ]*result[:Tlast])
             result[:HL]              = LOG2 / result[:Kel]
             result[:AUCinf]          = result[:AUClast] + result[:Clast] / result[:Kel]
+            result[:AUCinf_pred]     = result[:AUClast] + result[:Clast_pred] / result[:Kel]
             result[:AUCpct]          = (result[:AUCinf] - result[:AUClast]) / result[:AUCinf] * 100.0
         end
         #-----------------------------------------------------------------------
@@ -656,9 +671,21 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
             else
                 result[:Ctau] = data.obs[end]
             end
+            daucpartl  = 0.0
+            daumcpartl = 0.0
+
+            if data.time[ncas] != data.dosetime.time
+                result[:Ctaudose] = minimum(data.obs[ncas:ncae])
+                daucpartl, daumcpartl = aucpart(data.dosetime.time, data.time[ncas], result[:Ctaudose], data.obs[ncas], calcm, false)
+            end
+
+
+            #result[:Ctaudose] = minimum(data.obs[ncas+1:ncae])
+            #daucpartl, daumcpartl = aucpart(data.dosetime.time, data.time[ncas+1], result[:Ctaudose], data.obs[ncas+1], calcm, false)
+
             result[:Ctaumin]  = min(result[:Ctau], result[:Cdose], minimum(data.obs[ncas+1:ncae]))
-            result[:AUCtau]   = sum(aucpartl[pmask])  + eaucpartl
-            result[:AUMCtau]  = sum(aumcpartl[pmask]) + eaumcpartl
+            result[:AUCtau]   = sum(aucpartl[pmask])  + eaucpartl  + daucpartl
+            result[:AUMCtau]  = sum(aumcpartl[pmask]) + eaumcpartl + daumcpartl
             result[:Cavg]     = result[:AUCtau]/data.dosetime.tau
             if result[:Ctaumin] != 0
                 result[:Swing]    = (result[:Cmax] - result[:Ctaumin])/result[:Ctaumin]
@@ -822,30 +849,44 @@ Pharmacodynamics non-compartment analysis for PD subjects set.
 function nca!(data::DataSet{PDSubject}; verbose = false, io::IO = stdout)
         results = Array{PKPDProfile, 1}(undef, 0)
         for i = 1:length(data.data)
-            push!(results, nca!(data.data[i]))
+            push!(results, nca!(data.data[i]; verbose = verbose, io = io))
         end
         return DataSet(results)
 end
 
 
 """
-    nca!(data::PDSubject; verbose = false, io::IO = stdout)::PKPDProfile{PDSubject}
+    nca!(data::UPKSubject; verbose = false, io::IO = stdout)::PKPDProfile{UPKSubject}
 
 Pharmacodynamics non-compartment analysis for urine data.
 """
 function nca!(data::UPKSubject; verbose = false, io::IO = stdout)::PKPDProfile{UPKSubject}
 
     result   = Dict()
-    result[:midtime] = (data.stime .+ data.etime) / 2
-    result[:exrate]  = (data.obs .* data.vol) ./  (data.stime .+ data.etime)
-    result[:ar]      =  sum(data.obs .* data.vol)
-    result[:maxrate] = maximum(result[:exrate])
+    #time  = (data.stime .+ data.etime) / 2
+    #exrate   = (data.conc .* data.vol) ./  (data.etime .- data.stime)
+    result[:maxrate], result[:mTmax], result[:mTmaxn] = ctmax(data.time, data.obs)
+    result[:ar]      = sum(data.conc .* data.vol)
     result[:volume]  = sum(data.vol)
+
     #result[:AUCrate]
     #result[:tmaxrate]
     #result[:arp]
     #result[:Kel]
     #result[:HL]
+    return PKPDProfile(data, result; method = :upk)
+end
+"""
+    nca!(data::DataSet{UPKSubject}; verbose = false, io::IO = stdout)
+
+Pharmacodynamics non-compartment analysis for PD subjects set.
+"""
+function nca!(data::DataSet{UPKSubject}; verbose = false, io::IO = stdout)
+    results = Vector{PKPDProfile}(undef, 0)
+    for i = 1:length(data)
+        push!(results, nca!(data[i]; verbose = verbose, io = io))
+    end
+    return DataSet(results)
 end
 #=
 stime::Vector
@@ -919,9 +960,7 @@ function pkimport(data::DataFrame, sort::Array, rule::LimitRule; conc::Symbol, t
     end
     =#
     results = DataSet(results)
-    if rule.lloq !== NaN
-         applyncarule!(results, rule)
-    end
+    applyncarule!(results, rule)
     return results
 end
 """
@@ -1011,29 +1050,36 @@ function pdimport(data::DataFrame; resp::Symbol, time::Symbol, bl = 0, th = NaN)
     datai = sort(data[!,[time, resp]], time)
     return DataSet([PDSubject(datai[!, time], datai[!, resp], bl, th)])
 end
+#-------------------------------------------------------------------------------
+function upkimport(data::DataFrame, sort::Array; stime::Symbol, etime::Symbol, conc::Symbol, vol::Symbol)::DataSet
+    results  = Array{UPKSubject, 1}(undef, 0)
+    getdatai(data, sort, [stime, etime, conc, vol], (x, y) -> push!(results, UPKSubject(copy(x[!, stime]), copy(x[!, etime]), copy(x[!, conc]), copy(x[!, vol]); sort =Dict(sort .=> collect(y)))); sortby = stime)
+    results = DataSet(results)
+    return results
+end
 
 #-------------------------------------------------------------------------------
-    function applyncarule(data::PKSubject, rule::LimitRule)
+    function applyncarule!(data::PKSubject, rule::LimitRule)
         cmax, tmax, tmaxn = ctmax(data)
         #NaN Rule
         obsn = length(data)
         if rule.nan !== NaN
             for i = 1:length(data)
-                if data.obs[i] === NaN data.obs[i] = rule.nan end
+                if data.obs[i] === NaN || data.obs[i] === missing data.obs[i] = rule.nan end
             end
         end
-
         #LLOQ rule
-        for i = 1:obsn
-            if data.obs[i] < rule.lloq
-                if i <= tmaxn
-                    data.obs[i] = rule.btmax
-                else
-                    data.obs[i] = rule.atmax
+        if rule.lloq !== NaN
+            for i = 1:obsn
+                if data.obs[i] <= rule.lloq
+                    if i <= tmaxn
+                        data.obs[i] = rule.btmax
+                    else
+                        data.obs[i] = rule.atmax
+                    end
                 end
             end
         end
-
         #NaN Remove rule
         if rule.rm
             filterv   = data.obs .!== NaN
@@ -1041,18 +1087,19 @@ end
             data.obs  = data.obs[filterv]
         end
     end
+    #=
     function applyncarule(data::DataSet{PKSubject}, rule::LimitRule)
     end
-    #=
+
     function applyncarule(data::PKPDProfile{PKSubject}, rule::LimitRule)
     end
     function applyncarule(data::DataSet{PKPDProfile{PKSubject}}, rule::LimitRule)
     end
-    =#
 
-    function applyncarule!(data::PKSubject, rule::LimitRule)
+    function applyncarule(data::PKSubject, rule::LimitRule)
 
     end
+    =#
     function applyncarule!(data::DataSet{PKSubject}, rule::LimitRule)
         for i in data
             applyncarule!(i, rule)
diff --git a/src/plots.jl b/src/plots.jl
index 6db1ffa..8ad9f5c 100644
--- a/src/plots.jl
+++ b/src/plots.jl
@@ -47,7 +47,6 @@ end
 
     seriestype  --> :line
     xlabel      --> "Time"
-    ylabel      --> "Concentration"
     link        --> :both
     legend      --> true
     grid        --> false
@@ -73,7 +72,7 @@ function plotlabel(d)
     return title
 end
 
-function plot(subj::T; title = nothing, legend = true, label = "AUTO", xlims = nothing, plotstyle::PKPlotStyle = PKPLOTSTYLE[1]) where T <: AbstractSubject
+function plot(subj::T; title = nothing, legend = true, label = "AUTO", ylabel = "Concentration", xlims = nothing, plotstyle::PKPlotStyle = PKPLOTSTYLE[1]) where T <: AbstractSubject
 
     if title === nothing
         title = plotlabel(subj.sort)
@@ -87,6 +86,7 @@ function plot(subj::T; title = nothing, legend = true, label = "AUTO", xlims = n
         markercolor = plotstyle.markercolor,
         legend      = legend,
         label       = label,
+        ylabel      = ylabel,
         xlims       = xlims,
         )
     return p
@@ -115,7 +115,7 @@ function usort(data::DataSet{T}, list) where T <: AbstractData
     dl
 end
 
-function pageplot(pagedatatmp, styledict, typesort; title = "", legend = true, xlims = nothing)
+function pageplot(pagedatatmp, styledict, typesort; title = "", ylabel = "Concentration", legend = true, xlims = nothing)
     utypes    = keys(styledict)
     labels    = Vector{String}(undef, 0)
     fst       = true
@@ -135,7 +135,7 @@ function pageplot(pagedatatmp, styledict, typesort; title = "", legend = true, x
                     end
                 end
                 if fst
-                    p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = style, label = label)
+                    p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = style, label = label, ylabel = ylabel)
                     fst = false
                 else
                     plot!(                  d; legend = legend, xlims = xlims, plotstyle = style, label = label)
@@ -145,13 +145,13 @@ function pageplot(pagedatatmp, styledict, typesort; title = "", legend = true, x
     end
     p
 end
-function pageplot(pagedatatmp; title = "", legend = true, xlims = nothing)
+function pageplot(pagedatatmp; title = "", ylabel = "Concentration", legend = true, xlims = nothing)
     fst       = true
     p         = nothing
     for d in pagedatatmp
         label = "AUTO"
         if fst
-            p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = PKPLOTSTYLE[1], label = label)
+            p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = PKPLOTSTYLE[1], label = label, ylabel = ylabel)
             fst = false
         else
             plot!(                  d; legend = legend, xlims = xlims, plotstyle = PKPLOTSTYLE[1], label = label)
@@ -160,7 +160,7 @@ function pageplot(pagedatatmp; title = "", legend = true, xlims = nothing)
     p
 end
 
-function plot(data::DataSet{T}; title = nothing, legend = true, xlims = nothing, pagesort = nothing, typesort = nothing) where T <: AbstractSubject
+function plot(data::DataSet{T}; title = nothing, ylabel = "Concentration", legend = true, xlims = nothing, pagesort = nothing, typesort = nothing) where T <: AbstractSubject
     # Style types
 
     styledict  = nothing
@@ -199,17 +199,17 @@ function plot(data::DataSet{T}; title = nothing, legend = true, xlims = nothing,
                 end
             end
             if typesort !== nothing
-                p = pageplot(pagedatatmp, styledict, typesort; title = title, legend = legend, xlims = xlims)
+                p = pageplot(pagedatatmp, styledict, typesort; title = title, legend = legend, xlims = xlims, ylabel = ylabel)
             else
-                p = pageplot(pagedatatmp; title = title, legend = legend, xlims = xlims)
+                p = pageplot(pagedatatmp; title = title, legend = legend, xlims = xlims, ylabel = ylabel)
             end
             push!(plots, p)
         end
     else
         if typesort !== nothing
-            p = pageplot(data, styledict, typesort; title = title, legend = legend, xlims = xlims)
+            p = pageplot(data, styledict, typesort; title = title, legend = legend, xlims = xlims, ylabel = ylabel)
         else
-            p = pageplot(data; title = title, legend = legend, xlims = xlims)
+            p = pageplot(data; title = title, legend = legend, xlims = xlims, ylabel = ylabel)
         end
         push!(plots, p)
     end
diff --git a/test/csv/nanpk.csv b/test/csv/nanpk.csv
new file mode 100644
index 0000000..a60d43e
--- /dev/null
+++ b/test/csv/nanpk.csv
@@ -0,0 +1,22 @@
+Concentration,Time,Subject,Formulation
+,0,1,T
+0.2,1,1,T
+0.3,2,1,T
+0.4,3,1,T
+0.3,4,1,T
+,5,1,T
+,6,1,T
+0,0,2,T
+0.3,1,2,T
+0.4,2,2,T
+0.5,3,2,T
+,4,2,T
+0.1,5,2,T
+0.05,6,2,T
+0,0,1,R
+,1,1,R
+0.9,2,1,R
+0.8,3,1,R
+0.3,4,1,R
+0,5,1,R
+0.1,6,1,R
diff --git a/test/csv/upkdata.csv b/test/csv/upkdata.csv
new file mode 100644
index 0000000..636a2e9
--- /dev/null
+++ b/test/csv/upkdata.csv
@@ -0,0 +1,6 @@
+subj	conc	st	et	vol
+1	1	0	1	1
+1	2	1	2	2
+1	2	2	6	3
+1	1	6	12	3
+1	1	12	18	2
diff --git a/test/pktest.jl b/test/pktest.jl
index 072d75d..3eadd34 100644
--- a/test/pktest.jl
+++ b/test/pktest.jl
@@ -216,7 +216,7 @@ println(" ---------------------------------- ")
 
     #AUCtau
     #=
-    4947.052768
+    @test round.(df[!, :AUCtau], sigdigits = 6) == round.([4947.052768
     6265.637822
     2863.392745
     5008.644865
@@ -225,7 +225,7 @@ println(" ---------------------------------- ")
     4096.937951 (!)
     3802.30601
     4531.779637
-    3744.768624
+    3744.768624], sigdigits = 6)
     =#
 
     # TAU 0.25 - 9
@@ -245,9 +245,12 @@ println(" ---------------------------------- ")
     142.566125
     122.7865], sigdigits = 6)
 
-    #AUC Tau (!!!)
-    #=
-    1268.275631
+    pk   = ClinicalTrialUtilities.nca!(pkds, calcm = :luld, intp = :logt)
+    df   = DataFrame(pk; unst = true)
+    sort!(df, :Subject)
+
+    #AUC Tau
+    @test round.(df[!, :AUCtau], sigdigits = 6) == round.([1268.275631
     1831.820528
     754.6493604
     1336.480932
@@ -256,8 +259,8 @@ println(" ---------------------------------- ")
     1079.366854
     766.620245
     1219.631864
-    970.3062692
-    =#
+    970.3062692], sigdigits = 6)
+
 
     # TAU 0.0 - 100
     ClinicalTrialUtilities.setdosetime!(pkds, ClinicalTrialUtilities.DoseTime(dose = 100, time = 0.0, tau = 100))
@@ -277,7 +280,7 @@ println(" ---------------------------------- ")
     23.98401112], sigdigits = 6)
 
 #=
-    12646.63632
+    @test round.(df[!, :AUCtau], sigdigits = 6) == round.([12646.63632
     11996.6718
     7195.902904
     11794.11692
@@ -286,7 +289,7 @@ println(" ---------------------------------- ")
     8395.400098 (!)
     8930.999936 (!)
     10727.4135 (!)
-    6389.420453
+    6389.420453], sigdigits = 6)
 =#
 
     # Log-trapezoidal Linear Interpolation
@@ -359,8 +362,32 @@ println(" ---------------------------------- ")
     p    = ClinicalTrialUtilities.plot(pkds; pagesort = [:Date], typesort = [:Subject])
     @test length(p) == 2
     print(io, pk[1].subject.keldata)
+
+
+    # NaN PK LimitRule test
+    pkds = ClinicalTrialUtilities.pkimport(nanpkdata, [:Subject, :Formulation]; conc = :Concentration, time = :Time)
+    rule = ClinicalTrialUtilities.LimitRule(0, 0, NaN, 0, true)
+    ClinicalTrialUtilities.applyncarule!(pkds, rule)
+    @test length(pkds[1]) == 5
+    @test pkds[3].obs[1]  == 0.0
+    @test pkds[2].obs[5]  == 0.1
+    pkds = ClinicalTrialUtilities.pkimport(nanpkdata, [:Subject, :Formulation]; conc = :Concentration, time = :Time)
+    rule = ClinicalTrialUtilities.LimitRule(0, 0, NaN, -1, false)
+    ClinicalTrialUtilities.applyncarule!(pkds, rule)
+    @test pkds[3].obs[6]  === NaN
 end
 
+println(" ---------------------------------- ")
+@testset " UPK                    " begin
+    upk = ClinicalTrialUtilities.upkimport(upkdata, [:subj]; stime = :st, etime = :et, conc = :conc, vol = :vol)
+    unca = ClinicalTrialUtilities.nca!(upk[1])
+    unca = ClinicalTrialUtilities.nca!(upk)
+    @test unca[1][:maxrate] ≈ 4.0
+    @test unca[1][:mTmax]   ≈ 1.5
+    @test unca[1][:ar]      ≈ 16.0
+    @test unca[1][:volume]  ≈ 11.0
+     p    = ClinicalTrialUtilities.plot(upk, ylabel = "Excretion")
+end
 
 println(" ---------------------------------- ")
 @testset "  PD                    " begin
diff --git a/test/testdata.jl b/test/testdata.jl
index 22cf00a..2a1ffaa 100644
--- a/test/testdata.jl
+++ b/test/testdata.jl
@@ -1,5 +1,7 @@
 #Simple PK DataFrame
 pkdata         = CSV.File(path*"/csv/pkdata.csv") |> DataFrame
+#Simple PK DataFrame
+nanpkdata         = CSV.File(path*"/csv/nanpk.csv") |> DataFrame
 #Simple PD DataFrame
 pddata         = CSV.File(path*"/csv/pddata.csv") |> DataFrame
 # Multiple subjects PK DataFrame
@@ -8,6 +10,8 @@ pkdata2        = CSV.File(path*"/csv/pkdata2.csv") |> DataFrame
 #Pinheiro, J. C. and Bates, D. M. (2000), Mixed-Effects Models in S and S-PLUS, Springer, New York. (Appendix A.10)
 #Hand, D. and Crowder, M. (1996), Practical Longitudinal Data Analysis, Chapman and Hall, London.
 glucose2       = CSV.File(path*"/csv/glucose2.csv") |> DataFrame
+#Simple urine PK DataFrame
+upkdata         = CSV.File(path*"/csv/upkdata.csv") |> DataFrame
 #Simple frequency dataset
 freqdat        = CSV.File(path*"/csv/freqdat.csv") |> DataFrame
 #Small dataset with negative and zero value

From a1afa5b54438f9e6d546b58a512be60a01fada60 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Thu, 13 Feb 2020 18:08:25 +0300
Subject: [PATCH 10/15] PartAUC fix & validate

---
 src/pk.jl | 70 +++++++++++++++++++++++++++++++++----------------------
 1 file changed, 42 insertions(+), 28 deletions(-)

diff --git a/src/pk.jl b/src/pk.jl
index e87db60..601bf87 100644
--- a/src/pk.jl
+++ b/src/pk.jl
@@ -537,6 +537,8 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
 
         doseaucpart  = 0.0
         doseaumcpart = 0.0
+        daucpartl    = 0.0
+        daumcpartl   = 0.0
 
         for i = result[:Obsnum]:-1:1
             if data.obs[i] > 0
@@ -564,22 +566,18 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
         pmask   = Array{Bool, 1}(undef, result[:Obsnum] - 1)
         pmask  .= true
 
-        #sum all AUC parts
-        result[:AUCall]  = sum(aucpartl)
-        result[:AUMCall] = sum(aumcpartl)
-        #-----------------------------------------------------------------------
-        #-----------------------------------------------------------------------
-        #Exclude all AUC parts from observed concentation before 0 or less
-        #Need elaborate!!!!
-        for i = result[:Tmaxn]:result[:Obsnum] - 1
-            if data.obs[i+1] <= 0 * data.obs[i+1] pmask[i:end] .= false; break end
-        end
-        result[:AUClast]  = sum(aucpartl[pmask])
-        result[:AUMClast] = sum(aumcpartl[pmask])
-        result[:MRTlast]       = result[:AUMClast] / result[:AUClast]
+        ncas    = nothing
+        ncae    = nothing
 
-        #-----------------------------------------------------------------------
+        if  data.dosetime.time !== NaN && data.dosetime.tau > 0
+            ncas, ncae       = ncarange(data, data.dosetime.time, data.dosetime.tau)
+            result[:Ctaumin] = minimum(data.obs[ncas:ncae])
+            #if data.time[ncas] != data.dosetime.time
+            #    daucpartl, daumcpartl = aucpart(data.dosetime.time, data.time[ncas], result[:Ctaumin], data.obs[ncas], calcm, false)
+            #end
+        end
 
+        #sum all AUC parts
         if data.dosetime.time < data.time[1]
             #Dosetime < first time
         end
@@ -591,8 +589,13 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
                     break
                 else
                 #Cdose calculation
-                    if data.obs[i] > 0
-                        result[:Cdose] = linpredict(data.time[i] , data.time[i + 1], data.dosetime.time, data.obs[i], data.obs[i + 1])
+
+
+                    #if data.obs[i] > 0
+                        #result[:Cdose] = linpredict(data.time[i] , data.time[i + 1], data.dosetime.time, data.obs[i], data.obs[i + 1])
+                        #esult[:Cdose] = data.obs[i]
+                    if  data.dosetime.time !== NaN && data.dosetime.tau > 0
+                        result[:Cdose] = result[:Ctaumin]
                     else
                         result[:Cdose] = 0
                     end
@@ -605,6 +608,23 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
                 pmask[i]     = false
             end
         end
+
+        result[:AUCall]  = sum(aucpartl[pmask])  + doseaucpart
+        result[:AUMCall] = sum(aumcpartl[pmask]) + doseaumcpart
+        #-----------------------------------------------------------------------
+        #-----------------------------------------------------------------------
+        #Exclude all AUC parts from observed concentation before 0 or less
+        #Need elaborate!!!!
+        for i = result[:Tmaxn]:result[:Obsnum] - 1
+            if data.obs[i+1] <= 0 * data.obs[i+1] pmask[i:end] .= false; break end
+        end
+        result[:AUClast]  = sum(aucpartl[pmask])
+        result[:AUMClast] = sum(aumcpartl[pmask])
+        result[:MRTlast]       = result[:AUMClast] / result[:AUClast]
+
+        #-----------------------------------------------------------------------
+
+
         #-----------------------------------------------------------------------
         #Elimination
         keldata                = KelData()
@@ -651,7 +671,6 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
         #Steady-state
         tautime    = data.dosetime.time + data.dosetime.tau
         if data.dosetime.tau > 0
-            ncas, ncae = ncarange(data, data.dosetime.time, data.dosetime.tau)
             #tautime    = data.dosetime.time + data.dosetime.tau
             eaucpartl  = eaumcpartl = 0.0
             if tautime < data.time[end]
@@ -671,21 +690,16 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
             else
                 result[:Ctau] = data.obs[end]
             end
-            daucpartl  = 0.0
-            daumcpartl = 0.0
-
-            if data.time[ncas] != data.dosetime.time
-                result[:Ctaudose] = minimum(data.obs[ncas:ncae])
-                daucpartl, daumcpartl = aucpart(data.dosetime.time, data.time[ncas], result[:Ctaudose], data.obs[ncas], calcm, false)
-            end
+            #daucpartl  = 0.0
+            #daumcpartl = 0.0
 
 
             #result[:Ctaudose] = minimum(data.obs[ncas+1:ncae])
             #daucpartl, daumcpartl = aucpart(data.dosetime.time, data.time[ncas+1], result[:Ctaudose], data.obs[ncas+1], calcm, false)
 
-            result[:Ctaumin]  = min(result[:Ctau], result[:Cdose], minimum(data.obs[ncas+1:ncae]))
-            result[:AUCtau]   = sum(aucpartl[pmask])  + eaucpartl  + daucpartl
-            result[:AUMCtau]  = sum(aumcpartl[pmask]) + eaumcpartl + daumcpartl
+            #esult[:Ctaumin]  = min(result[:Ctau], result[:Cdose], minimum(data.obs[ncas+1:ncae]))
+            result[:AUCtau]   = sum(aucpartl[pmask])  + eaucpartl  + doseaucpart
+            result[:AUMCtau]  = sum(aumcpartl[pmask]) + eaumcpartl + doseaumcpart
             result[:Cavg]     = result[:AUCtau]/data.dosetime.tau
             if result[:Ctaumin] != 0
                 result[:Swing]    = (result[:Cmax] - result[:Ctaumin])/result[:Ctaumin]
@@ -699,7 +713,7 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io:
             end
             result[:Fluc]     = (result[:Cmax] - result[:Ctaumin])/result[:Cavg]
             result[:Fluctau]  = (result[:Cmax] - result[:Ctau])/result[:Cavg]
-            #result[:AccInd]
+            result[:Accind]   = 1 / (1 - (exp(-result[:Kel] * data.dosetime.tau)))
         end
         if verbose
             aucpartlsum  = similar(aucpartl)

From 238629fe1308d0a6d42a4ce1ce4975e15471b583 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Thu, 13 Feb 2020 18:08:33 +0300
Subject: [PATCH 11/15] update

---
 cange.log      |   5 +++
 test/pktest.jl | 104 +++++++++++++++++++++++++++++++++++++++++++------
 2 files changed, 98 insertions(+), 11 deletions(-)

diff --git a/cange.log b/cange.log
index 8cd7bf9..0b3fb7f 100644
--- a/cange.log
+++ b/cange.log
@@ -1,6 +1,11 @@
 v0.2.3
 
   - ConTab struct
+  - Urine PK
+  - NCA Tau PK partial areas fix
+  - NCA test update
+  - Plotting improve
+  - minor fix
 
 
 v0.2.2
diff --git a/test/pktest.jl b/test/pktest.jl
index 3eadd34..12f087b 100644
--- a/test/pktest.jl
+++ b/test/pktest.jl
@@ -199,7 +199,7 @@ println(" ---------------------------------- ")
     # Linear Up Log Down - Log Interpolation
     # TAU 0 - 36
     ClinicalTrialUtilities.setdosetime!(pkds, ClinicalTrialUtilities.DoseTime(dose = 100, time = 0, tau = 36))
-    pk   = ClinicalTrialUtilities.nca!(pkds, calcm = :luldt, intp = :logt)
+    pk   = ClinicalTrialUtilities.nca!(pkds, calcm = :luld, intp = :logt)
     df   = DataFrame(pk; unst = true)
     sort!(df, :Subject)
 
@@ -214,23 +214,35 @@ println(" ---------------------------------- ")
     113.3749669
     72.57153458], sigdigits = 6)
 
+    @test round.(df[!, :AUClast], sigdigits = 6) == round.([9573.810559
+    10054.28648
+    5392.457219
+    9297.096334
+    9519.180874
+    6948.985621
+    6988.772632
+    7073.092214
+    8303.358586
+    5486.838889], sigdigits = 6)
+
+
     #AUCtau
-    #=
+
     @test round.(df[!, :AUCtau], sigdigits = 6) == round.([4947.052768
     6265.637822
     2863.392745
     5008.644865
     5415.246398
     3882.203934
-    4096.937951 (!)
+    4096.937951
     3802.30601
     4531.779637
     3744.768624], sigdigits = 6)
-    =#
+
 
     # TAU 0.25 - 9
     ClinicalTrialUtilities.setdosetime!(pkds, ClinicalTrialUtilities.DoseTime(dose = 100, time = 0.25, tau = 9))
-    pk   = ClinicalTrialUtilities.nca!(pkds, calcm = :luldt, intp = :logt)
+    pk   = ClinicalTrialUtilities.nca!(pkds, calcm = :luld, intp = :logt)
     df   = DataFrame(pk; unst = true)
     sort!(df, :Subject)
 
@@ -262,9 +274,67 @@ println(" ---------------------------------- ")
     970.3062692], sigdigits = 6)
 
 
+    @test round.(df[!, :AUCall], sigdigits = 6) == round.([9566.596809
+    10054.28648
+    5392.457219
+    9297.096334
+    9519.180874
+    6948.985621
+    6988.772632
+    7058.818964
+    8302.368086
+    5486.838889], sigdigits = 6)
+    #Swing tau
+    @test round.(df[!, :Swingtau], sigdigits = 6) == round.([0.38590184
+    0.27426721
+    0.240584404
+    0.468637128
+    0.057960085
+    0.445494022
+    0.240196282
+    0.628152237
+    0.173820219
+    0.021952739], sigdigits = 6)
+
+    #Fluctuation%_Tau
+
+    @test round.(df[!, :Fluctau] * 100, sigdigits = 6) == round.([37.71452462
+    27.61899665
+    24.36421266
+    44.8121175
+    6.369060133
+    38.4975876
+    24.74919661
+    62.65259796
+    18.28649133
+    2.500189968], sigdigits = 6)
+
+    @test round.(df[!, :Kel], sigdigits = 6) == round.([0.003384744
+    0.014106315
+    0.00329143
+    0.007695344
+    0.006813328
+    0.007692281
+    0.012458956
+    0.00893008
+    0.005645865
+    0.017189737], sigdigits = 6)
+
+    @test round.(df[!, :Accind], sigdigits = 6) == round.([33.3295745
+    8.38726982
+    34.26016611
+    14.94451581
+    16.81301567
+    14.95026394
+    9.427514161
+    12.94903929
+    20.18432092
+    6.976692409], sigdigits = 6)
+
+
     # TAU 0.0 - 100
     ClinicalTrialUtilities.setdosetime!(pkds, ClinicalTrialUtilities.DoseTime(dose = 100, time = 0.0, tau = 100))
-    pk   = ClinicalTrialUtilities.nca!(pkds, calcm = :luldt, intp = :logt)
+    pk   = ClinicalTrialUtilities.nca!(pkds, calcm = :luld, intp = :logt)
     df   = DataFrame(pk; unst = true)
     sort!(df, :Subject)
 
@@ -279,18 +349,29 @@ println(" ---------------------------------- ")
     80.05171762
     23.98401112], sigdigits = 6)
 
-#=
+    @test round.(df[!, :AUClast], sigdigits = 6) == round.([9573.810559
+    10054.28648
+    5392.457219
+    9297.096334
+    9519.180874
+    6948.985621
+    6988.772632
+    7073.092214
+    8303.358586
+    5486.838889], sigdigits = 6)
+
+
     @test round.(df[!, :AUCtau], sigdigits = 6) == round.([12646.63632
     11996.6718
     7195.902904
     11794.11692
     12274.83395
     8729.151856
-    8395.400098 (!)
-    8930.999936 (!)
-    10727.4135 (!)
+    8395.400098
+    8930.999936
+    10727.4135
     6389.420453], sigdigits = 6)
-=#
+
 
     # Log-trapezoidal Linear Interpolation
 
@@ -321,6 +402,7 @@ println(" ---------------------------------- ")
     24849.129
     7940.0834], sigdigits = 6)
 
+    #LULDT
     pk   = ClinicalTrialUtilities.nca!(pkds, calcm = :luldt, io = io, verbose = true)
     df   = DataFrame(pk; unst = true)
     sort!(df, :Subject)

From a982bf40bca31d1b6ac35487db9c3b21b1edd5a7 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Thu, 13 Feb 2020 18:10:30 +0300
Subject: [PATCH 12/15] update

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 6a18237..f0301aa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 authors = ["Vladimir Arnautov (mail@pharmcat.net)"]
 name = "ClinicalTrialUtilities"
 uuid = "535c2557-d7d0-564d-8ff9-4ae146c18cfe"
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

From 0a58fb52ab1372f770301e842252d9b69ec718eb Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Thu, 13 Feb 2020 18:12:24 +0300
Subject: [PATCH 13/15] fix

---
 src/freque.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/freque.jl b/src/freque.jl
index d13e642..44ca667 100644
--- a/src/freque.jl
+++ b/src/freque.jl
@@ -82,3 +82,8 @@ end
 function fisher(t::ConTab{2, 2})
     fisher(t.tab)
 end
+
+#=
+function StatsBase.confint(t::ConTab{2, 2})
+end
+=#

From 642fa61d326234f3158be27de14377fabf266281 Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Thu, 13 Feb 2020 18:29:32 +0300
Subject: [PATCH 14/15] fix

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index f0301aa..1c196ce 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,6 +30,6 @@ DataFrames = "0.19, 0.20"
 QuadGK = "2.0, 2.1"
 SpecialFunctions = "0.8, 0.9"
 Roots = "0.7, 0.8"
-CSV = "0.5"
+CSV = "0.4, 0.5"
 RecipesBase = "0.7"
 Reexport = "0.1, 0.2"

From aacf162f6deaa7d0af33939560726a3e8333b8cd Mon Sep 17 00:00:00 2001
From: PharmCat <v.s.arnautov@yandex.ru>
Date: Thu, 13 Feb 2020 18:34:22 +0300
Subject: [PATCH 15/15] fix

---
 Project.toml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 1c196ce..f73e521 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,6 +13,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 
 [extras]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -30,6 +31,6 @@ DataFrames = "0.19, 0.20"
 QuadGK = "2.0, 2.1"
 SpecialFunctions = "0.8, 0.9"
 Roots = "0.7, 0.8"
-CSV = "0.4, 0.5"
+CSV = "0.5"
 RecipesBase = "0.7"
 Reexport = "0.1, 0.2"