SSD cannot handle arbitrary contact IDs
fhagemann opened this issue · 3 comments
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
.
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
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:
- Solve the initial issue by making sure that we
- either force the contacts to be labelled from
1
toN
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.) - 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 theWeightingPotentials
:
via the index in theArray
(or whatever we might store them in) or via the contact IDs.
I agree that #306 is overkill.
- either force the contacts to be labelled from
- 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.
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).