check_convergence.jl

Convergence is checked starting from the simulated third cardiac cycle. Waveforms are first loaded and subdivided depending on the number of cycles already passed. This task is accomplished by function extractWaveform. Waveforms are then resampled in order to contain a smaller number of points. This would speed up the following error calculation. These two steps are taken by resampleAndCalculateError function. Loading, re-sampling, and error calculation operations are done for all the quantities calculated by openBF. All these functions are wrapped into checkAllQuantities which is called inside main.jl loop.

function checkAllQuantities \(\rightarrow\) max_error::Float

Parameters:
v ::Vessel vessel of which the error is being checked.
passed_cycles ::Int number of cardiac cycles already passed.
n_pts ::Int number of points in the resampled waveform. 

function checkAllQuantities(v             :: Vessel,
                            passed_cycles :: Int64,
                            n_pts         :: Int64)
Functioning
er is an array which will contain the percentage error for each quantity in the current vessel. qs collection contains all the suffixes needed to build .out filenames for each quantity. 
  er = zeros(Float64, 5)

  qs = ["_P.out", "_Q.out", "_A.out", "_c.out", "_u.out"]

For each quantity q listed in qs, the .out filename is created by joining the label inside v and the suffix q. Then the desired waveforms are loaded by extractWaveform, re-sampled, and the error is calculated by resampleAndCalculateError. Eventually, the maximum error scored among all the quantities is returned.

Returns:
max_error ::Float maximum error within all the quantities in the current vessel. 
  idx = 1
  for q in qs
    filename = join([v.label, q])

    W = extractWaveform(filename, passed_cycles)

    er[idx] = resampleAndCalculateError(v, W, passed_cycles, n_pts)

    idx += 1
  end

  return maximum(er)
end

function checkAllQuantities \(\rightarrow\) max_error::Float

Parameters:
vs ::Array{Vessel, 1} collection of vessels.
passed_cycles ::Int number of cardiac cycles already passed.
n_pts ::Int number of points in the resampled waveform.
Functioning
To compute the error for all the vessels, checkAllQuantities is re-defined to handle also a vector of Vessels. checkAllQuantities is called recursively.
Returns:
max_error ::Float maximum error among all the vessels. 

function checkAllQuantities(vs            :: Array{Vessel, 1},
                            passed_cycles :: Int64,
                            n_pts         :: Int64)

  ers = zeros(Float64, length(vs))

  idx = 1
  for v in vs

    ers[idx] = checkAllQuantities(v, passed_cycles, n_pts)
    idx += 1

  end

  return maximum(ers)
end

function extractWaveform \(\rightarrow\) Fx::Array{Array{Float, 2}}

Parameters:
filename ::String current time in seconds.
cycles ::Int \(\Delta t\), current time step. 

function extractWaveform(filename :: String, cycles :: Int64)
Functioning
filename is a string referring to the .out file containing the waveform to be imported. This is read by readdlm function and stored into F array. F is a timesteps \(\times\) 6 array; the first column contains the time variable, the remaining columns are one for each monitored node along the current vessel. 
  F = readdlm(filename)

These three lines find the cardiac period t by knowing how many cycles have already passed.

 
  sf1 = int(length(F[:,1])*(cycles-1)/cycles)
  T1 = F[sf1:end,1] .- F[sf1,1]
  t = T1[end]

Fx is an empty array which will contain one array for each cardiac cycle.

 
  Fx = Array{Float64, 2}[]

This loop skims through F, separates waveforms, and add them to Fx. si is the waveform source index, the index from which the waveform starts. ti is the waveform terminal index, the index at which the waveforms ends because a new cardiac cycle is due to start.

 
  si = 1
  ti = 1
  idx = 1
  for i in 1:length(F[:,1])-1
    if F[i,1] < t*idx && F[i+1, 1] >= t*idx
      ti = i
      push!(Fx, F[si:ti,:])
      si = ti
      idx += 1
    end
  end
Returns:
Fx ::Array{Array{Float, 2}} 
  return Fx
end

function resampleAndCalculateError \(\rightarrow\) max_error::Float

Parameters:
W ::Array{Array{Float, 2}} collection of all waveforms during the passed cardiac cycles for the current quantity.
passed_cycles ::Int number of cardiac cycles already simulated.
n_pts ::Int number of points to be contained in the re-sampled waveform. 

function resampleAndCalculateError(v, W :: Array{Array{Float64, 2}, 1},
                                   passed_cycles :: Int64,
                                   n_pts         :: Int64)
Functioning
F will contain the re-sampled waveform, er is an empty variable for the maximum error. 
  F = zeros(Float64, n_pts, 2)
  er = 0

The error is computed for all the monitor nodes. Their values are stored in W’s columns from 2 to 6. Each waveform is re-sampled by using resample function from scipy.signal library. In order to compute the error, only the last and the previous waveforms are analysed.

 
  for j in 2:6

    F[:,1], t = signal.resample(W[end-1][:,j], n_pts, W[end-1][:,1])
    F[:,2], t = signal.resample(W[end  ][:,j], n_pts, W[end  ][:,1])

    er = abs(F[:,1] .- F[:,2])./F[:,1]*100
  end
Returns:
max_error ::Float maximum error on the current quantity along the the current vessel. 
  return maximum(er)
end

function checkConvergence(vessels :: Array{Vessel, 1},
  convergence_tollerance :: Float64)

  for v in vessels
    err = checkConvergence(v, convergence_tollerance)
    if err > convergence_tollerance
      return err
    end
  end
  return err
end

function checkConvergence(v :: Vessel,
  convergence_tollerance :: Float64)

  qs = ["_P", "_Q"]

  for q in qs
    filename_last = join([v.label, q, ".last"])
    filename_temp = join([v.label, q, ".temp"])

    err = checkQuantity(filename_last, filename_temp, convergence_tollerance)
    if err > convergence_tollerance
      return err
    end
  end
  return err
end

function checkQuantity(last, temp, convergence_tollerance)

  s1 = readdlm(last)
  s2 = readdlm(temp)
  n_pts = 1000

  r1 = n_pts/length(s1[:,2])
  r2 = n_pts/length(s2[:,2])
  for j in 2:5
    sr1 = resample(s1[:,j], r1)
    sr2 = resample(s2[:,j], r2)
    for i in 20:length(sr1)-20
      err = abs(100*(sr2[i]-sr1[i])/sr1[i])
      if err > convergence_tollerance
        return err
      end
    end
  end
  return 9.99
end