SSD cannot handle arbitrary contact IDs

[fhagemann]

Right now, a config file can only be read in, if the contacts of the detector are enumerated in order from 1 to N
(where N is the number of contacts).

Choosing something different from this results in an error when trying to display the Simulation:

detectors:
  - semiconductor:
      material: HPGe
      geometry: 
        # ...
    contacts:
      - material: HPGe
        potential: 2000
        id: 101
        geometry:
          # ...
      - material: HPGe
        potential: 0
        id: 102
        geometry:
          # ...

sim = Simulation("config_file.yaml")

Simulation{Float32, Cartesian} - Coordinate system: Cartesian
  Environment Material: Vacuum
  Detector: ExampleCuboid
  Electric potential: missing
  Charge density: missing
  Impurity scale: missing
  Fix Charge density: missing
  Dielectric distribution: missing
  Point types: missing
  Electric field: missing
  Weighting potentials: 
    Contact 101: Error showing value of type Simulation{Float32, Cartesian}:
ERROR: BoundsError: attempt to access 2-element Vector{Any} at index [101]
Stacktrace:
  [1] getindex(A::Vector{Any}, i1::Int64)
    @ Base ./array.jl:861
  [2] println(io::IOContext{Base.TTY}, sim::Simulation{Float32, Cartesian})
    @ SolidStateDetectors ~/.julia/packages/SolidStateDetectors/ETZbh/src/Simulation/Simulation.jl:134
  [3] show
    @ ~/.julia/packages/SolidStateDetectors/ETZbh/src/Simulation/Simulation.jl:142 [inlined]
  [4] show(io::IOContext{Base.TTY}, #unused#::MIME{Symbol("text/plain")}, sim::Simulation{Float32, Cartesian})
    @ SolidStateDetectors ~/.julia/packages/SolidStateDetectors/ETZbh/src/Simulation/Simulation.jl:145
  [5] (::REPL.var"#43#44"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:266
  [6] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:510
  [7] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:259
  [8] display(d::REPL.REPLDisplay, x::Any)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:271
  [9] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:328
 [10] #invokelatest#2
    @ ./essentials.jl:716 [inlined]
 [11] invokelatest
    @ ./essentials.jl:714 [inlined]
 [12] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:293
 [13] (::REPL.var"#45#46"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:277
 [14] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:510
 [15] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:275
 [16] (::REPL.var"#do_respond#66"{Bool, Bool, REPL.var"#77#87"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:846
 [17] #invokelatest#2
    @ ./essentials.jl:716 [inlined]
 [18] invokelatest
    @ ./essentials.jl:714 [inlined]
 [19] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/LineEdit.jl:2493
 [20] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL /opt/julia-1.7.2/share/julia/stdlib/v1.7/REPL/src/REPL.jl:1232
 [21] (::REPL.var"#49#54"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ./task.jl:423

We also get an error running the simulation:

simulate!(sim)

ERROR: BoundsError: attempt to access 2-element Vector{Any} at index [101]
Stacktrace:
 [1] setindex!
   @ ./essentials.jl:479 [inlined]
 [2] apply_initial_state!(sim::Simulation{Float32, Cartesian}, ::Type{WeightingPotential}, contact_id::Int64, grid::Grid{Float32, 3, Cartesian, Tuple{SolidStateDetectors.DiscreteAxis{Float32, :infinite, :infinite, IntervalSets.ClosedInterval{Float32}}, SolidStateDetectors.DiscreteAxis{Float32, :infinite, :infinite, IntervalSets.ClosedInterval{Float32}}, SolidStateDetectors.DiscreteAxis{Float32, :infinite, :infinite, IntervalSets.ClosedInterval{Float32}}}}; not_only_paint_contacts::Bool, paint_contacts::Bool, depletion_handling::Bool)
   @ SolidStateDetectors ~/.julia/packages/SolidStateDetectors/ETZbh/src/Simulation/Simulation.jl:539
 [3] _calculate_potential!(sim::Simulation{Float32, Cartesian}, potential_type::UnionAll, contact_id::Int64; grid::Missing, convergence_limit::Float64, refinement_limits::Vector{Float64}, min_tick_distance::Missing, max_tick_distance::Missing, max_distance_ratio::Int64, depletion_handling::Bool, use_nthreads::Int64, sor_consts::Missing, max_n_iterations::Int64, n_iterations_between_checks::Int64, not_only_paint_contacts::Bool, paint_contacts::Bool, verbose::Bool, device_array_type::Type{Array})
   @ SolidStateDetectors ~/.julia/packages/SolidStateDetectors/ETZbh/src/Simulation/Simulation.jl:932
 [4] #calculate_weighting_potential!#228
   @ ~/.julia/packages/SolidStateDetectors/ETZbh/src/Simulation/Simulation.jl:1070 [inlined]
 [5] simulate!(sim::Simulation{Float32, Cartesian}; convergence_limit::Float64, refinement_limits::Vector{Float64}, min_tick_distance::Missing, max_tick_distance::Missing, max_distance_ratio::Int64, depletion_handling::Bool, use_nthreads::Int64, sor_consts::Missing, max_n_iterations::Int64, device_array_type::Type{Array}, not_only_paint_contacts::Bool, paint_contacts::Bool, verbose::Bool)
   @ SolidStateDetectors ~/.julia/packages/SolidStateDetectors/ETZbh/src/Simulation/Simulation.jl:1299
 [6] simulate!(sim::Simulation{Float32, Cartesian})
   @ SolidStateDetectors ~/.julia/packages/SolidStateDetectors/ETZbh/src/Simulation/Simulation.jl:1283
 [7] top-level scope
   @ REPL[6]:1

It might be great to allow for arbitrary contact IDs (I guess that was implemented in one of the very first versions of SSD),
especially for users with a Python background who would like to give their first contact the ID 0 instead of 1.

[lmh91]

So we have currently two PRs, #290 and #306 facing this issue differently.

Regarding the contact_id issue in the Event struct. I think we don't have to deal with that.
I would like to remove the Event struct anyhow soon and replace it by an Table
which we partially already have: see simulate_waveforms.
It will just be an MCTable with just one unique event id in the event_id column.
The channel id column will hold the respect contact id.

So we just have to deal with how to store the weighting potentials.
I think creating an custom array type like in #306 is a bit overkill.
The array is also not type stable as the types of the weighting potentials might be different.
It is also a lot of code which I think is not needed. Its basically an own implementation of a dictionary.

So I would favor the Dict approach here even though I am also not a fan.
We also just could use a Table and do it via NamedTuples:

sim.weighting_potentials = Table(
    contact_id = [1, 200, 34],
    weighting_potential = [wp_1, wp_200, wp_34]
)

Another possibility (which could go along with one of the other approaches) is to define

struct ContactID
     id::Int
end

Such that we can overload certain getindex methods.
However, then one would always have to wrap the id into that type:
sim.weighting_potentials[ContactID(1)]

Please share opinions on this @oschulz @fhagemann @schustermartin @SebastianRuffert

[fhagemann]

So, the main issue right now is that SSD cannot handle config files in which the contacts are not labelled from 1 to N according to their appearance in the config file.

In addition, the current implementation is not type stable (it is just a Union of Missing or Array{WeightingPotential}).

I think we need to address these two things separately:

1. Solve the initial issue by making sure that we
   1. either force the contacts to be labelled from 1 to N and throw an error or warning if this is not the case
      (this might not be favoured as we might wanna combine config files via splitting.)
   2. or allow for arbitrary IDs and adapt the code where needed to make this possible.
      If we choose this option, we should clarify how the user accesses the WeightingPotentials:
      via the index in the Array (or whatever we might store them in) or via the contact IDs.
      I agree that Introduce CustomIDArray to access Arrays with contact IDs #306 is overkill.
2. Think of a type stable solution on how to store WeightingPotentials.

In my opinion, we should realize 1.i in a small PR (to remove the bug from the code for now) and think about a long-term solution that covers 1.ii and 2. I am fine with using both Dict or Table for storing the WeightingPotentials.

sim.weighting_potentials = Table(
    contact_id = [1, 200, 34],
    weighting_potential = [wp_1, wp_200, wp_34]
)

However, this does also not really seem type stable to me.

[lmh91]

I don't think type stability actually matters here.
The mutable struct Simulation is only supposed to be a very high level collection object to keep everything together for the user. And we won't use sim::Simulation in performance relevant code anyhow. For example:

wp = sim.weighting_potentials[1] # <- This is fine if it is not type stable.
for event in evt_table
    get_signal(event, wp) 
end

The only way (as far as I know) of making it entirely really type stabile would be to use a Tuple of all weighting potentials.
However, this would create a lot of code generation and would be really bad if there would be some highly pixelated detector (e.g. more than 50 contacts).

I can imagine some use cases where it might makes sense for the user to specify contact IDs which are not labelled from 1 to N.

I agree with your opinion with realizing 1.i for now and then think about a long-term solution a bit more.
Maybe in combination with the refactoring of the Event interface (#226).