From 04b36875c9e2db0679bc367793254a728595b366 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Tue, 3 Oct 2023 18:53:39 -0400 Subject: [PATCH 01/18] mc gen model --- .gitignore | 1 + Manifest.toml | 9 +++++---- Project.toml | 1 + src/GalileoEvents.jl | 1 + src/gms/gms.jl | 4 +++- src/gms/mc_gm.jl | 18 +++++++++++------- src/utils/distributions.jl | 8 ++++---- src/utils/ramp.obj | 15 +++++++++++++++ src/utils/scenes.jl | 20 ++++++++++++-------- test/gms/mc_gm.jl | 36 ++++++++++++++++++++++++------------ 10 files changed, 77 insertions(+), 36 deletions(-) create mode 100644 src/utils/ramp.obj diff --git a/.gitignore b/.gitignore index 634cdc0..92f5024 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,4 @@ worker-* tmp* dask-worker-space __pycache__ +env.d/jenv/ \ No newline at end of file diff --git a/Manifest.toml b/Manifest.toml index 22b471e..5317cf9 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,13 +2,13 @@ julia_version = "1.9.3" manifest_format = "2.0" -project_hash = "8ead7856fbb2572c3c0a941b94dca0886303c698" +project_hash = "3fee985c9eca9c1165906a81b33e7ff0c8e89188" [[deps.Accessors]] -deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Requires", "Test"] -git-tree-sha1 = "954634616d5846d8e216df1298be2298d55280b2" +deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Test"] +git-tree-sha1 = "a7055b939deae2455aa8a67491e034f735dd08d3" uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" -version = "0.1.32" +version = "0.1.33" [deps.Accessors.extensions] AccessorsAxisKeysExt = "AxisKeys" @@ -19,6 +19,7 @@ version = "0.1.32" [deps.Accessors.weakdeps] AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5" IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" + Requires = "ae029012-a4dd-5104-9daa-d747884805df" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" diff --git a/Project.toml b/Project.toml index 79b2d73..2a186be 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["belledon"] version = "0.1.0" [deps] +Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" diff --git a/src/GalileoEvents.jl b/src/GalileoEvents.jl index 77b4588..8edea2b 100644 --- a/src/GalileoEvents.jl +++ b/src/GalileoEvents.jl @@ -1,5 +1,6 @@ module GalileoEvents +using Accessors using Gen using Gen_Compose using PhySMC diff --git a/src/gms/gms.jl b/src/gms/gms.jl index 34e24fe..db12e91 100644 --- a/src/gms/gms.jl +++ b/src/gms/gms.jl @@ -1,3 +1,5 @@ +using FillArrays + export GMParams, GMState, Material, @@ -56,7 +58,7 @@ $(TYPEDSIGNATURES) A uniform prior over given materials """ -function MaterialPrior(ms::Vector{Material}) +function MaterialPrior(ms::Vector{<: Material}) n = length(ms) ws = Fill(1.0 / n, n) MaterialPrior(ms, ws) diff --git a/src/gms/mc_gm.jl b/src/gms/mc_gm.jl index d139915..dae9ea2 100644 --- a/src/gms/mc_gm.jl +++ b/src/gms/mc_gm.jl @@ -17,11 +17,12 @@ $(TYPEDFIELDS) struct MCParams <: GMParams # prior material_prior::MaterialPrior - physics_prior::Vector{PhysPrior} + physics_prior::PhysPrior # simulation sim::BulletSim template::BulletState n_objects::Int64 + obs_noise::Float64 end """ @@ -29,8 +30,9 @@ $(TYPEDSIGNATURES) Initializes `MCParams` from a constructed scene in pybullet. """ -function MCParams(client::PyObject, objs::Vector{PyObject}, - mprior::MaterialPrior, pprior::Vector{PhysPrior}) +function MCParams(client::Int64, objs::Vector{Int64}, + mprior::MaterialPrior, pprior::PhysPrior, + obs_noise::Float64=0.) # configure simulator with the provided # client id sim = BulletSim(;client=client) @@ -41,7 +43,7 @@ function MCParams(client::PyObject, objs::Vector{PyObject}, # Note: alternative latents will be suggested by the `prior` template = BulletState(sim, rigid_bodies) - MCParams(mprior, pprior, sim, template, length(objs)) + MCParams(mprior, pprior, sim, template, length(objs), obs_noise) end struct MCState <: GMState @@ -56,9 +58,9 @@ end @gen function mc_object_prior(ls::RigidBodyLatents, gm::MCParams) # sample material - mi = @trace(categorical(gm.material_weights), :material) + mi = @trace(categorical(gm.material_prior.material_weights), :material) # sample physical properties - phys_params = gm.phys_params[mi] + phys_params = gm.physics_prior mass_mu, mass_sd = phys_params.mass mass = @trace(trunc_norm(mass_mu, mass_sd, 0., Inf), :mass) @@ -94,7 +96,8 @@ end @gen function kernel(t::Int, prev_state::MCState, gm::MCParams) sim_step::BulletState = PhySMC.step(gm.sim, prev_state.bullet_state) - obs = @trace(Gen.Map(observe_position)(sim_step.kinematics), :observe) + noises = Fill(gm.obs_noise, length(sim_step.kinematics)) + obs = @trace(Gen.Map(observe_position)(sim_step.kinematics, noises), :observe) next_state = MCState(sim_step) return next_state end @@ -102,6 +105,7 @@ end @gen function mc_gm(t::Int, gm::MCParams) init_state = @trace(mc_prior(gm), :prior) # simulate `t` timesteps + println(init_state) states = @trace(Gen.Unfold(kernel)(t, init_state, gm), :kernel) return states end diff --git a/src/utils/distributions.jl b/src/utils/distributions.jl index c8dc7b2..c597f82 100644 --- a/src/utils/distributions.jl +++ b/src/utils/distributions.jl @@ -8,14 +8,14 @@ struct NoisyMatrix <: Gen.Distribution{Array{Float64}} end const mat_noise = NoisyMatrix() -function Gen.logpdf(::NoisyMatrix, x::Array{Float64}, mu::Array{U}, noise::T) where {U<:Real,T<:Real} +function Gen.logpdf(::NoisyMatrix, x::Array{Float64}, mu::Array{<:Real}, noise::T) where {T<:Real} var = noise * noise diff = x - mu vec = diff[:] return -(vec' * vec)/ (2.0 * var) - 0.5 * log(2.0 * pi * var) end; -function Gen.random(::NoisyMatrix, mu::Array{U}, noise::T) where {U<:Real,T<:Real} +function Gen.random(::NoisyMatrix, mu::Array{<:Real}, noise::T) where {T<:Real} mat = copy(mu) for i in CartesianIndices(mu) mat[i] = mu[i] + randn() * noise @@ -30,7 +30,7 @@ struct LogUniform <: Gen.Distribution{Float64} end const log_uniform = LogUniform() -function Gen.logpdf(::LogUniform, x::Float64, low::T, high::T) where {U<:Real,T<:Real} +function Gen.logpdf(::LogUniform, x::Float64, low::T, high::T) where {T<:Real} l = log(low) h = log(high) v = log(x) @@ -38,7 +38,7 @@ function Gen.logpdf(::LogUniform, x::Float64, low::T, high::T) where {U<:Real,T< return (v >= l && v <= h) ? -log(h-l) : -Inf end -function Gen.random(::LogUniform, low::T, high::T) where {U<:Real,T<:Real} +function Gen.random(::LogUniform, low::T, high::T) where {T<:Real} d = uniform(log(low), log(high)) exp(d) end diff --git a/src/utils/ramp.obj b/src/utils/ramp.obj new file mode 100644 index 0000000..59fd2ad --- /dev/null +++ b/src/utils/ramp.obj @@ -0,0 +1,15 @@ +v 0.000000 0.000000 0.000000 + +v 1.000000 0.000000 0.000000 + +v 1.000000 1.000000 0.000000 + +v 0.000000 1.000000 0.000000 + +v 0.000000 0.000000 1.000000 + +v 0.000000 1.000000 1.000000 + +f 1 2 3 4 + +f 1 2 5 6 \ No newline at end of file diff --git a/src/utils/scenes.jl b/src/utils/scenes.jl index 954d1c1..fc1f3dc 100644 --- a/src/utils/scenes.jl +++ b/src/utils/scenes.jl @@ -6,8 +6,11 @@ export ramp obj_positions::NTuple{2}; slope, ramp_intersection) """ function ramp( + mass_ratio::Float64, + obj_frictions::NTuple{2, Float64} = (.5, .5), + obj_positions::NTuple{2, Float64} = (0.5, 1.5), slope::Float64=2/3, - tableRampIntersection::Float64=0., + tableRampIntersection::Float64=0. ) # for debugging #client = @pycall pb.connect(pb.GUI)::Int64 @@ -56,7 +59,8 @@ function ramp( end # add a ramp - ramp_col_id = pb.createCollisionShape(pb.GEOM_MESH, fileName="examples/ramp/ramp.obj", physicsClientId=client, meshScale=[2, base_dims[2], slope*2]) + pb.setAdditionalSearchPath("/project") + ramp_col_id = pb.createCollisionShape(pb.GEOM_MESH, fileName="src/utils/ramp.obj", physicsClientId=client, meshScale=[2, base_dims[2], slope*2]) ramp_position = [-2+tableRampIntersection, -base_dims[2]/2, 0] ramp_obj_id = pb.createMultiBody(baseCollisionShapeIndex=ramp_col_id, basePosition=ramp_position, physicsClientId=client) pb.changeDynamics(ramp_obj_id, -1; mass=0.0, restitution=0.9, physicsClientId=client) @@ -90,18 +94,18 @@ function ramp( obj_on_ramp_col_id = pb.createCollisionShape(pb.GEOM_BOX, halfExtents=obj_ramp_dims/2, physicsClientId=client) lift = obj_ramp_dims[3]/2 position = [ - -1+tableRampIntersection+lift*cos(theta_radians), + -2+2*obj_positions[1]+tableRampIntersection+lift*cos(theta_radians), 0, - 1*slope-lift*sin(theta_radians) + (2-2*obj_positions[1])*slope-lift*sin(theta_radians) ] obj_on_ramp_obj_id = pb.createMultiBody(baseCollisionShapeIndex=obj_on_ramp_col_id, basePosition=position, baseOrientation=orientation, physicsClientId=client) - pb.changeDynamics(obj_on_ramp_obj_id, -1; mass=1.0, restitution=0.9, physicsClientId=client) - + pb.changeDynamics(obj_on_ramp_obj_id, -1; mass=mass_ratio, restitution=0.9, lateralFriction=obj_frictions[1], physicsClientId=client) + # add an object on the table that will collide with the object on the ramp as that one slides down obj_on_table_dims = [0.2, 0.2, 0.1] obj_on_table_col_id = pb.createCollisionShape(pb.GEOM_BOX, halfExtents=obj_on_table_dims/2, physicsClientId=client) - obj_on_table_obj_id = pb.createMultiBody(baseCollisionShapeIndex=obj_on_table_col_id, basePosition=[1, 0, obj_on_table_dims[3]/2], physicsClientId=client) - pb.changeDynamics(obj_on_table_obj_id, -1; mass=1.0, restitution=0.9, physicsClientId=client) + obj_on_table_obj_id = pb.createMultiBody(baseCollisionShapeIndex=obj_on_table_col_id, basePosition=[2.5*(obj_positions[2]-1), 0, obj_on_table_dims[3]/2], physicsClientId=client) + pb.changeDynamics(obj_on_table_obj_id, -1; mass=1.0, restitution=0.9, lateralFriction=obj_frictions[2], physicsClientId=client) (client, obj_on_ramp_obj_id, obj_on_table_obj_id) end diff --git a/test/gms/mc_gm.jl b/test/gms/mc_gm.jl index d00e39a..2bb0374 100644 --- a/test/gms/mc_gm.jl +++ b/test/gms/mc_gm.jl @@ -3,21 +3,20 @@ using GalileoEvents mass_ratio = 2.0 obj_frictions = (0.3, 0.3) -obj_positions = () # TODO fill in... +obj_positions = (0.5, 1.5) mprior = MaterialPrior([unknown_material]) pprior = PhysPrior((3.0, 10.0), # mass (0.5, 10.0), # friction (0.2, 1.0)) # restitution +obs_noise = 0.05 +t = 120 -t = 60 - -# TODO: this should evaluate without errors function forward_test() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) - mc_params = MCParams(client, [a,b], mprior, pprior) - trace, _ = Gen.generate(t, mc_gm) - display(get_choices(trace)) + mc_params = MCParams(client, [a,b], mprior, pprior, obs_noise) + trace, _ = Gen.generate(mc_gm, (t, mc_params)) + #display(get_choices(trace)) end @@ -25,13 +24,26 @@ end # compare the final positions (they should be different) function update_test() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) - mc_params = MCParams(client, [a,b], mprior, pprior) - trace, _ = Gen.generate(t, mc_gm) + mc_params = MCParams(client, [a,b], mprior, pprior, obs_noise) + trace, _ = Gen.generate(mc_gm, (t, mc_params)) addr = :prior => :objects => 1 => :mass - cm = Gen.choicemap(addr => trace[addr] + 3.0) - trace2 = Gen.update(trace, cm) - + cm = Gen.choicemap(addr => trace[addr] + 3) + trace2, _ = Gen.update(trace, cm) # compare final positions + for i in 1:2 + @assert get_submap(get_choices(trace), :kernel=>120=>:observe=>i) != get_submap(get_choices(trace2), :kernel=>120=>:observe=>i) + end + + return trace, trace2 end + +function main() + forward_test() + update_test() +end + +if abspath(PROGRAM_FILE) == @__FILE__ + main() +end \ No newline at end of file From 735af0cf796666dec2d9b16ee7fa33d298442945 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 4 Oct 2023 11:48:23 -0400 Subject: [PATCH 02/18] compare retval for update --- test/gms/mc_gm.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/gms/mc_gm.jl b/test/gms/mc_gm.jl index 2bb0374..5257958 100644 --- a/test/gms/mc_gm.jl +++ b/test/gms/mc_gm.jl @@ -20,8 +20,6 @@ function forward_test() end -# TODO: use `Gen.update` to change a traces physical latents and -# compare the final positions (they should be different) function update_test() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) mc_params = MCParams(client, [a,b], mprior, pprior, obs_noise) @@ -32,9 +30,10 @@ function update_test() trace2, _ = Gen.update(trace, cm) # compare final positions - for i in 1:2 - @assert get_submap(get_choices(trace), :kernel=>120=>:observe=>i) != get_submap(get_choices(trace2), :kernel=>120=>:observe=>i) - end + t=120 + pos1 = Vector(get_retval(trace)[t].bullet_state.kinematics[1].position) + pos2 = Vector(get_retval(trace2)[t].bullet_state.kinematics[1].position) + @assert pos1 != pos2 return trace, trace2 end From 5ca5c37886700f5a6db6b6b5e30cbab100e963cc Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 4 Oct 2023 14:11:06 -0400 Subject: [PATCH 03/18] change point model WIP --- src/gms/cp_gm_pb.jl | 171 ++++++++++++++++++++++++++++++++++++++++++++ src/gms/mc_gm.jl | 1 - 2 files changed, 171 insertions(+), 1 deletion(-) create mode 100644 src/gms/cp_gm_pb.jl diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl new file mode 100644 index 0000000..f649099 --- /dev/null +++ b/src/gms/cp_gm_pb.jl @@ -0,0 +1,171 @@ + +export CPParams, + CPState, + CPModel + +using LinearAlgebra:norm + +## Generative Model + components + +struct ObjectBelief + congruent::Bool + physical_latents::RigidBodyLatents +end + +abstract type EventRelation end + +struct Collision <: EventRelation + a::RigidBody + b::RigidBody +end + +struct EventState + active_events::Vector{EventRelation} +end + +""" +Parameters for change point model + +$(TYPEDEF) + +--- + +$(TYPEDFIELDS) +""" +struct CPParams <: GMParams + # prior + material_prior::MaterialPrior + physics_prior::PhysPrior + # event relations + event_concepts::Vector{Type{EventRelation}} + # simulation + sim::BulletSim + template::BulletState + n_objects::Int64 + obs_noise::Float64 +end + +struct CPState <: GMState + bullet_state::BulletState + event_state::EventState +end + +## PRIOR + +@gen function cp_object_prior(ls::RigidBodyLatents, gm::CPParams) + # sample material + mi = @trace(categorical(gm.material_prior.material_weights), :material) + + # sample physical properties + phys_params = gm.physics_prior + mass_mu, mass_sd = phys_params.mass + mass = @trace(trunc_norm(mass_mu, mass_sd, 0., Inf), :mass) + fric_mu, fric_sd = phys_params.friction + friction = @trace(trunc_norm(fric_mu,fric_sd, 0., 1.), :friction) + res_low, res_high = phys_params.restitution + restitution = @trace(uniform(res_low, res_high), :restitution) + + # package + new_ls = setproperties(ls.data; + mass = mass, + lateralFriction = friction, + restitution = restitution) + new_latents::RigidBodyLatents = RigidBodyLatents(new_ls) + return new_latents +end + +@gen function cp_prior(params::CPParams) + # initialize the kinematic state + latents = params.template.latents + params_filled = Fill(params, length(latents)) + new_latents = @trace(Gen.Map(cp_object_prior)(latents, params_filled), :objects) + bullet_state = setproperties(params.template; latents = new_latents) + + # initialize the event state + event_state = EventState([]) + + init_state = CPState(bullet_state, event_state) + return init_state +end + +""" +Bernoulli weight that event relation holds +""" +function predicate(Type{Collision}, x::RigidBodyState, y::RigidBodyState) + d = norm(x.position, y.position) # l2 distance + clamp(1.0 - d, 0., 1.) +end + +# in case of collision: Gaussian drift update of mass and restitution +@gen static function _collision_clause(event_idx, latents::RigidBodyLatents) + new_mass = @trace(trunc_norm(latents.mass, .1, 0., Inf), :mass) + new_restitution = @trace(trunc_norm(latents.restitution, .1, res_low, res_high), :restitution) + return set_properties(latents; mass = new_mass, restitution = new_restitution) +end + +""" +Returns the generative function over latents for the event relation +""" +function clause(Type{Collision}) + _collision_clause +end + + +function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelation}}) + return event_concepts + for EventRelation in event_concepts + # TODO: decide if valid + end +end + +@gen function update_latents(state, latents) + new_latents = setproperties(latents.data; + mass = mass, + lateralFriction = friction, + restitution = restitution) + new_latents::RigidBodyLatents = RigidBodyLatents(new_latents) + return new_latents +end + +@gen function event_kernel(state::CPState, event_concepts::Vector{Type{EventRelation}}) + # filter out invalid event relations (e.g., a collision is already active between a and b) + relations = valid_relations(state, event_concepts) + + # pairwise weights of event concepts + weights = pairwise(event_concepts, state) + event_idx = @trace(categorical(weights), :event_idx) + + new_latents = @trace(Switch(map(clause, event_concepts))(event_idx, state.bullet_state), :event) + update_latents(state, new_latents) + return active_events +end + +# for one object, observe the noisy positions +@gen function observe_position(k::RigidBodyState, noise::Float64) + pos = k.position + # add noise to position + obs = @trace(broadcasted_normal(pos, noise), :positions) + return obs +end + +# for one time step, run event and physics kernel +@gen (static) function kernel(t::Int, prev_state::CPState, params::CPParams) + # event kernel + event_state = @trace(event_kernel(prev_state, params.event_concepts), :events) + + # simulate physics for the next step based on the current state and observe positions + bullet_state::BulletState = PhySMC.step(params.sim, state.bullet_state) + obs = @trace(Gen.Map(observe_position)(bullet_state.kinematics, noises), :observe) + + next_state = CPState(bullet_state, event_state) + return next_state +end + +@gen (static) function CPModel(t::Int, params::CPParams) + # initalize the kinematic and event state + init_state = @trace(cp_prior(params), :prior) + + # unfold the event and kinematic state over time + states = @trace(Gen.Unfold(kernel)(t, init_state, params), :kernel) + return states +end diff --git a/src/gms/mc_gm.jl b/src/gms/mc_gm.jl index dae9ea2..ab38f66 100644 --- a/src/gms/mc_gm.jl +++ b/src/gms/mc_gm.jl @@ -105,7 +105,6 @@ end @gen function mc_gm(t::Int, gm::MCParams) init_state = @trace(mc_prior(gm), :prior) # simulate `t` timesteps - println(init_state) states = @trace(Gen.Unfold(kernel)(t, init_state, gm), :kernel) return states end From e49f58ec521a7b6d88d36e7e327952fd48bf681e Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 4 Oct 2023 15:26:19 -0400 Subject: [PATCH 04/18] agenda for cp model --- src/gms/cp_gm_pb.jl | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index f649099..613af68 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -1,7 +1,7 @@ export CPParams, CPState, - CPModel + cp_model using LinearAlgebra:norm @@ -19,9 +19,6 @@ struct Collision <: EventRelation b::RigidBody end -struct EventState - active_events::Vector{EventRelation} -end """ Parameters for change point model @@ -47,7 +44,7 @@ end struct CPState <: GMState bullet_state::BulletState - event_state::EventState + active_events::Vector{EventRelation} end ## PRIOR @@ -113,6 +110,7 @@ end function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelation}}) return event_concepts + # TODO: replace by map in the end for EventRelation in event_concepts # TODO: decide if valid end @@ -127,16 +125,24 @@ end return new_latents end +# iterate over event concepts and evaluate predicates to active/deactive @gen function event_kernel(state::CPState, event_concepts::Vector{Type{EventRelation}}) # filter out invalid event relations (e.g., a collision is already active between a and b) - relations = valid_relations(state, event_concepts) + #relations = valid_relations(state, event_concepts) - # pairwise weights of event concepts - weights = pairwise(event_concepts, state) + # activate new events + # map active events to weights 3D tensor for birth decision + # evaluate predicates object-pairwise + weights = #call predicates event_concepts for each event_idx = @trace(categorical(weights), :event_idx) - new_latents = @trace(Switch(map(clause, event_concepts))(event_idx, state.bullet_state), :event) + # Switch combinator to evaluate clauses for each event, currently only one + new_latents = @trace(Gen.Switch(map(clause, event_concepts))(event_idx, state.bullet_state), :event) + # alternative: update latents with Map update_latents(state, new_latents) + + # some new events kill old events + return active_events end @@ -161,7 +167,7 @@ end return next_state end -@gen (static) function CPModel(t::Int, params::CPParams) +@gen (static) function cp_model(t::Int, params::CPParams) # initalize the kinematic and event state init_state = @trace(cp_prior(params), :prior) From c15287feb7be1017059abc5efd3c502eb7eca6c0 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 4 Oct 2023 15:40:22 -0400 Subject: [PATCH 05/18] add structure to changepoint model --- src/gms/cp_gm_pb.jl | 53 +++++++++++++++++++++++++++++---------------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 613af68..7397d57 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -5,13 +5,11 @@ export CPParams, using LinearAlgebra:norm -## Generative Model + components - -struct ObjectBelief - congruent::Bool - physical_latents::RigidBodyLatents -end +## Changepoint Model + components +""" +Event types +""" abstract type EventRelation end struct Collision <: EventRelation @@ -22,12 +20,6 @@ end """ Parameters for change point model - -$(TYPEDEF) - ---- - -$(TYPEDFIELDS) """ struct CPParams <: GMParams # prior @@ -42,6 +34,9 @@ struct CPParams <: GMParams obs_noise::Float64 end +""" +Current state of the change point model, simulation state and event state +""" struct CPState <: GMState bullet_state::BulletState active_events::Vector{EventRelation} @@ -49,6 +44,9 @@ end ## PRIOR +""" +initalizes prior beliefs about mass, friction and resitution of the given objects +""" @gen function cp_object_prior(ls::RigidBodyLatents, gm::CPParams) # sample material mi = @trace(categorical(gm.material_prior.material_weights), :material) @@ -71,6 +69,9 @@ end return new_latents end +""" +initializes belief about all objects and events +""" @gen function cp_prior(params::CPParams) # initialize the kinematic state latents = params.template.latents @@ -79,9 +80,9 @@ end bullet_state = setproperties(params.template; latents = new_latents) # initialize the event state - event_state = EventState([]) + active_events = Vector{EventRelation}() - init_state = CPState(bullet_state, event_state) + init_state = CPState(bullet_state, active_events) return init_state end @@ -93,7 +94,9 @@ function predicate(Type{Collision}, x::RigidBodyState, y::RigidBodyState) clamp(1.0 - d, 0., 1.) end -# in case of collision: Gaussian drift update of mass and restitution +""" +in case of collision: Gaussian drift update of mass and restitution# +""" @gen static function _collision_clause(event_idx, latents::RigidBodyLatents) new_mass = @trace(trunc_norm(latents.mass, .1, 0., Inf), :mass) new_restitution = @trace(trunc_norm(latents.restitution, .1, res_low, res_high), :restitution) @@ -107,7 +110,9 @@ function clause(Type{Collision}) _collision_clause end - +""" +objects that are already involved in some events should not be involved in new events +""" function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelation}}) return event_concepts # TODO: replace by map in the end @@ -116,6 +121,7 @@ function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelati end end + @gen function update_latents(state, latents) new_latents = setproperties(latents.data; mass = mass, @@ -125,7 +131,9 @@ end return new_latents end -# iterate over event concepts and evaluate predicates to active/deactive +""" +iterate over event concepts and evaluate predicates to active/deactive +""" @gen function event_kernel(state::CPState, event_concepts::Vector{Type{EventRelation}}) # filter out invalid event relations (e.g., a collision is already active between a and b) #relations = valid_relations(state, event_concepts) @@ -146,7 +154,9 @@ end return active_events end -# for one object, observe the noisy positions +""" +for one object, observe the noisy position in every dimension +""" @gen function observe_position(k::RigidBodyState, noise::Float64) pos = k.position # add noise to position @@ -154,7 +164,9 @@ end return obs end -# for one time step, run event and physics kernel +""" +for one time step, run event and physics kernel +""" @gen (static) function kernel(t::Int, prev_state::CPState, params::CPParams) # event kernel event_state = @trace(event_kernel(prev_state, params.event_concepts), :events) @@ -167,6 +179,9 @@ end return next_state end +""" +generate physical scene with changepoints in the belief state +""" @gen (static) function cp_model(t::Int, params::CPParams) # initalize the kinematic and event state init_state = @trace(cp_prior(params), :prior) From 8ce2d2a8aed838d9f79f2210d421dc546db0a18b Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Tue, 10 Oct 2023 22:04:23 -0400 Subject: [PATCH 06/18] continue changepoint model WIP --- Manifest.toml | 7 ++- Project.toml | 2 + src/gms/cp_gm_pb.jl | 122 ++++++++++++++++++++++++++++++-------------- src/gms/gms.jl | 2 +- test/gms/cp_gm.jl | 45 ++++++++++++++++ 5 files changed, 139 insertions(+), 39 deletions(-) create mode 100644 test/gms/cp_gm.jl diff --git a/Manifest.toml b/Manifest.toml index 5317cf9..598bd57 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.9.3" manifest_format = "2.0" -project_hash = "3fee985c9eca9c1165906a81b33e7ff0c8e89188" +project_hash = "6c92244b6b2c576cb0a17ab8a9dd1f5a2f86d403" [[deps.Accessors]] deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Test"] @@ -102,6 +102,11 @@ git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" version = "0.12.10" +[[deps.Combinatorics]] +git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860" +uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +version = "1.0.2" + [[deps.CommonSubexpressions]] deps = ["MacroTools", "Test"] git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" diff --git a/Project.toml b/Project.toml index 2a186be..0ce58f8 100644 --- a/Project.toml +++ b/Project.toml @@ -5,11 +5,13 @@ version = "0.1.0" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" +Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" Gen = "ea4f424c-a589-11e8-07c0-fd5c91b9da4a" Gen_Compose = "c1ef4dca-b0a6-4a35-b24b-46cbf3979a16" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" PhyBullet = "63daae69-5b14-439d-ac6f-096429ca839b" PhySMC = "79c1e2f5-7911-41a0-b248-4858717ddd79" diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 7397d57..525d0f9 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -1,9 +1,12 @@ export CPParams, CPState, - cp_model + cp_model, + EventRelation, + Collision using LinearAlgebra:norm +using Combinatorics ## Changepoint Model + components @@ -17,6 +20,8 @@ struct Collision <: EventRelation b::RigidBody end +struct NoEvent <: EventRelation end + """ Parameters for change point model @@ -26,7 +31,7 @@ struct CPParams <: GMParams material_prior::MaterialPrior physics_prior::PhysPrior # event relations - event_concepts::Vector{Type{EventRelation}} + event_concepts::Vector{Type{<:EventRelation}} # simulation sim::BulletSim template::BulletState @@ -34,12 +39,29 @@ struct CPParams <: GMParams obs_noise::Float64 end +function CPParams(client::Int64, objs::Vector{Int64}, + mprior::MaterialPrior, pprior::PhysPrior, + event_concepts::Vector{Type{<:EventRelation}}, + obs_noise::Float64=0.) + # configure simulator with the provided + # client id + sim = BulletSim(;client=client) + # These are the objects of interest in the scene + rigid_bodies = RigidBody.(objs) + # Retrieve the default latents for the objects + # as well as their initial positions + # Note: alternative latents will be suggested by the `prior` + template = BulletState(sim, rigid_bodies) + + CPParams(mprior, pprior, event_concepts, sim, template, length(objs), obs_noise) +end + """ Current state of the change point model, simulation state and event state """ struct CPState <: GMState bullet_state::BulletState - active_events::Vector{EventRelation} + active_events::Vector{Int64} end ## PRIOR @@ -72,7 +94,7 @@ end """ initializes belief about all objects and events """ -@gen function cp_prior(params::CPParams) +@gen (static) function cp_prior(params::CPParams) # initialize the kinematic state latents = params.template.latents params_filled = Fill(params, length(latents)) @@ -89,27 +111,51 @@ end """ Bernoulli weight that event relation holds """ -function predicate(Type{Collision}, x::RigidBodyState, y::RigidBodyState) - d = norm(x.position, y.position) # l2 distance +function predicate(t::Type{Collision}, a::RigidBodyState, b::RigidBodyState) + d = norm(Vector(a.position)-Vector(b.position)) # l2 distance clamp(1.0 - d, 0., 1.) end +""" +update latents of a single element +""" +@gen function update_latents(latents::BulletElemLatents) + new_mass = @trace(trunc_norm(latents.data.mass, .1, 0., Inf), :mass) + new_restitution = @trace(trunc_norm(latents.data.restitution, .1, 0., 1.), :restitution) + + new_latents = setproperties(latents.data; + mass = new_mass, + restitution = new_restitution) + new_latents = RigidBodyLatents(new_latents) + return new_latents +end + + + """ in case of collision: Gaussian drift update of mass and restitution# """ -@gen static function _collision_clause(event_idx, latents::RigidBodyLatents) - new_mass = @trace(trunc_norm(latents.mass, .1, 0., Inf), :mass) - new_restitution = @trace(trunc_norm(latents.restitution, .1, res_low, res_high), :restitution) - return set_properties(latents; mass = new_mass, restitution = new_restitution) +@gen (static) function _collision_clause(pair_idx::Vector{Int64}, latents::Vector{<:BulletElemLatents}) + new_latents = @trace(update_latents(latents[pair_idx[1]]), :new_latents_a) + new_latents = @trace(update_latents(latents[pair_idx[2]]), :new_latents_b) + return Type{Collision} +end + +@gen (static) function _no_event_clause(pair_idx, latents::Vector{BulletElemLatents}) + return Type{NoEvent} end """ Returns the generative function over latents for the event relation """ -function clause(Type{Collision}) +function clause(::Type{Collision}) _collision_clause end +function clause(::Type{NoEvent}) + _no_event_clause +end + """ objects that are already involved in some events should not be involved in new events """ @@ -121,35 +167,37 @@ function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelati end end - -@gen function update_latents(state, latents) - new_latents = setproperties(latents.data; - mass = mass, - lateralFriction = friction, - restitution = restitution) - new_latents::RigidBodyLatents = RigidBodyLatents(new_latents) - return new_latents -end - """ iterate over event concepts and evaluate predicates to active/deactive """ -@gen function event_kernel(state::CPState, event_concepts::Vector{Type{EventRelation}}) - # filter out invalid event relations (e.g., a collision is already active between a and b) +@gen function event_kernel(state::CPState, event_concepts::Vector{Type{<:EventRelation}}) + # TODO: filter out invalid event relations (e.g., a collision is already active between a and b) #relations = valid_relations(state, event_concepts) # activate new events - # map active events to weights 3D tensor for birth decision - # evaluate predicates object-pairwise - weights = #call predicates event_concepts for each + object_pairs = collect(combinations(state.bullet_state.kinematics, 2)) # get pairs of objects + # map active events to weights 2D tensor for birth decision + + weights = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs] + push!(weights, 1 - min(1, sum(weights))) # for no event + weights = weights ./ sum(weights) + #println(weights) + + pair_idx = repeat(collect(combinations(1:length(state.bullet_state.kinematics), 2)), length(event_concepts)) + push!(pair_idx, Vector([0,0])) # for no event + #println(pair_idx) + + # draw one event overall event_idx = @trace(categorical(weights), :event_idx) + active_events = [state.active_events..., event_idx] + #push!(active_events, event_idx) + #println(active_events) - # Switch combinator to evaluate clauses for each event, currently only one - new_latents = @trace(Gen.Switch(map(clause, event_concepts))(event_idx, state.bullet_state), :event) - # alternative: update latents with Map - update_latents(state, new_latents) + # Switch combinator to evaluate clauses for each event + events = vcat([repeat([event_type], length(object_pairs)) for event_type in event_concepts]..., NoEvent) + event = @trace(Gen.Switch(map(clause, events)...)(event_idx, pair_idx[event_idx], state.bullet_state.latents), :event) - # some new events kill old events + # TODO: some new events kill old events return active_events end @@ -157,7 +205,7 @@ end """ for one object, observe the noisy position in every dimension """ -@gen function observe_position(k::RigidBodyState, noise::Float64) +@gen (static) function observe_position(k::RigidBodyState, noise::Float64) pos = k.position # add noise to position obs = @trace(broadcasted_normal(pos, noise), :positions) @@ -169,13 +217,13 @@ for one time step, run event and physics kernel """ @gen (static) function kernel(t::Int, prev_state::CPState, params::CPParams) # event kernel - event_state = @trace(event_kernel(prev_state, params.event_concepts), :events) + active_events = @trace(event_kernel(prev_state, params.event_concepts), :events) # simulate physics for the next step based on the current state and observe positions - bullet_state::BulletState = PhySMC.step(params.sim, state.bullet_state) - obs = @trace(Gen.Map(observe_position)(bullet_state.kinematics, noises), :observe) + bullet_state::BulletState = PhySMC.step(params.sim, prev_state.bullet_state) + obs = @trace(Gen.Map(observe_position)(bullet_state.kinematics, Fill(params.obs_noise, params.n_objects)), :observe) - next_state = CPState(bullet_state, event_state) + next_state = CPState(bullet_state, active_events) return next_state end @@ -185,7 +233,7 @@ generate physical scene with changepoints in the belief state @gen (static) function cp_model(t::Int, params::CPParams) # initalize the kinematic and event state init_state = @trace(cp_prior(params), :prior) - + # unfold the event and kinematic state over time states = @trace(Gen.Unfold(kernel)(t, init_state, params), :kernel) return states diff --git a/src/gms/gms.jl b/src/gms/gms.jl index db12e91..164d556 100644 --- a/src/gms/gms.jl +++ b/src/gms/gms.jl @@ -80,4 +80,4 @@ struct PhysPrior end include("mc_gm.jl") -# include("cp_gm.jl") +include("cp_gm_pb.jl") diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl new file mode 100644 index 0000000..53080d7 --- /dev/null +++ b/test/gms/cp_gm.jl @@ -0,0 +1,45 @@ +using Gen +using GalileoEvents + +mass_ratio = 2.0 +obj_frictions = (0.3, 0.3) +obj_positions = (0.5, 1.5) + +mprior = MaterialPrior([unknown_material]) +pprior = PhysPrior((3.0, 10.0), # mass + (0.5, 10.0), # friction + (0.2, 1.0)) # restitution + +obs_noise = 0.05 +t = 120 + +function forward_test() + client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) + event_concepts = Type{<:EventRelation}[Collision] + cp_params = CPParams(client, [a,b, b], mprior, pprior, event_concepts, obs_noise) + trace, _ = Gen.generate(cp_model, (t, cp_params)) + display(get_choices(trace)) +end + + +function update_test() + client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) + mc_params = MCParams(client, [a,b], mprior, pprior, obs_noise) + trace, _ = Gen.generate(mc_gm, (t, mc_params)) + + addr = :prior => :objects => 1 => :mass + cm = Gen.choicemap(addr => trace[addr] + 3) + trace2, _ = Gen.update(trace, cm) + + # compare final positions + t=120 + pos1 = Vector(get_retval(trace)[t].bullet_state.kinematics[1].position) + pos2 = Vector(get_retval(trace2)[t].bullet_state.kinematics[1].position) + @assert pos1 != pos2 + + return trace, trace2 +end + + +forward_test() +#update_test() \ No newline at end of file From 082e2ab2a686d256eef1a5918378feaf6e5140b9 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 18 Oct 2023 16:11:23 -0400 Subject: [PATCH 07/18] experiment with changepoint model --- Manifest.toml | 2 +- src/gms/cp_gm_pb.jl | 132 ++++++++++++++++++++++++++++++-------------- test/gms/cp_gm.jl | 95 +++++++++++++++++++++++++++++-- 3 files changed, 180 insertions(+), 49 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 598bd57..7bcc38b 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -362,7 +362,7 @@ version = "0.72.9+1" [[deps.Gen]] deps = ["Compat", "DataStructures", "Distributions", "ForwardDiff", "FunctionalCollections", "JSON", "LinearAlgebra", "MacroTools", "Parameters", "Random", "ReverseDiff", "SpecialFunctions"] -git-tree-sha1 = "9878ff4ab1990f5647e89b4228a3c9da5f0e69c7" +path = "/home/dg963/GalileoEvents/env.d/jenv/dev/Gen" uuid = "ea4f424c-a589-11e8-07c0-fd5c91b9da4a" version = "0.4.6" diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 525d0f9..f673c07 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -61,7 +61,7 @@ Current state of the change point model, simulation state and event state """ struct CPState <: GMState bullet_state::BulletState - active_events::Vector{Int64} + active_events::Set{Int64} end ## PRIOR @@ -94,7 +94,7 @@ end """ initializes belief about all objects and events """ -@gen (static) function cp_prior(params::CPParams) +@gen function cp_prior(params::CPParams) # initialize the kinematic state latents = params.template.latents params_filled = Fill(params, length(latents)) @@ -102,7 +102,7 @@ initializes belief about all objects and events bullet_state = setproperties(params.template; latents = new_latents) # initialize the event state - active_events = Vector{EventRelation}() + active_events = Set{EventRelation}() init_state = CPState(bullet_state, active_events) return init_state @@ -113,9 +113,14 @@ Bernoulli weight that event relation holds """ function predicate(t::Type{Collision}, a::RigidBodyState, b::RigidBodyState) d = norm(Vector(a.position)-Vector(b.position)) # l2 distance - clamp(1.0 - d, 0., 1.) + @show d + clamp(1.0 - exp(-d), 0., 1.) end +# TODO: surface distances, use pybullet (contact maybe not helpful) + +# gen functional collections + """ update latents of a single element """ @@ -130,19 +135,17 @@ update latents of a single element return new_latents end - - """ in case of collision: Gaussian drift update of mass and restitution# """ -@gen (static) function _collision_clause(pair_idx::Vector{Int64}, latents::Vector{<:BulletElemLatents}) - new_latents = @trace(update_latents(latents[pair_idx[1]]), :new_latents_a) - new_latents = @trace(update_latents(latents[pair_idx[2]]), :new_latents_b) - return Type{Collision} +@gen function _collision_clause(pair_idx::Vector{Int64}, latents::Vector{BulletElemLatents}) + latents[pair_idx[1]] = @trace(update_latents(latents[pair_idx[1]]), :new_latents_a) + latents[pair_idx[2]] = @trace(update_latents(latents[pair_idx[2]]), :new_latents_b) + return latents end -@gen (static) function _no_event_clause(pair_idx, latents::Vector{BulletElemLatents}) - return Type{NoEvent} +@gen function _no_event_clause(pair_idx, latents::Vector{BulletElemLatents}) + return latents end """ @@ -157,7 +160,7 @@ function clause(::Type{NoEvent}) end """ -objects that are already involved in some events should not be involved in new events + """ function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelation}}) return event_concepts @@ -167,45 +170,90 @@ function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelati end end +## TODO: add Pkg revise + + +# map mcmc kernel (link: https://www.gen.dev/docs/stable/ref/mcmc/#Composing-Stationary-Kernels) +# set of proposal functions + +# change randomly clause choices of mass +# revise for which event, modify the categorical +# e.g. another event type for the same pair or another event for one event type + + +# cognition: computation around events +# modelling: changepoint model with Switch combinator +# inference: change types of events to help inference being faster + + """ iterate over event concepts and evaluate predicates to active/deactive """ @gen function event_kernel(state::CPState, event_concepts::Vector{Type{<:EventRelation}}) - # TODO: filter out invalid event relations (e.g., a collision is already active between a and b) - #relations = valid_relations(state, event_concepts) - - # activate new events - object_pairs = collect(combinations(state.bullet_state.kinematics, 2)) # get pairs of objects - # map active events to weights 2D tensor for birth decision - - weights = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs] - push!(weights, 1 - min(1, sum(weights))) # for no event - weights = weights ./ sum(weights) - #println(weights) - + println("-------") + object_pairs = collect(combinations(state.bullet_state.kinematics, 2)) pair_idx = repeat(collect(combinations(1:length(state.bullet_state.kinematics), 2)), length(event_concepts)) - push!(pair_idx, Vector([0,0])) # for no event + pair_idx = [[0,0], pair_idx...] # for no event #println(pair_idx) - - # draw one event overall - event_idx = @trace(categorical(weights), :event_idx) - active_events = [state.active_events..., event_idx] - #push!(active_events, event_idx) - #println(active_events) - # Switch combinator to evaluate clauses for each event - events = vcat([repeat([event_type], length(object_pairs)) for event_type in event_concepts]..., NoEvent) - event = @trace(Gen.Switch(map(clause, events)...)(event_idx, pair_idx[event_idx], state.bullet_state.latents), :event) + # map possible events to weight vector for birth decision using the predicates + predicates = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs] + print("Predicates: ") + println(predicates) - # TODO: some new events kill old events + # randomly draw a starting and an ending event + # birth + weights = copy(predicates) + for idx in state.active_events # active events should not be born again + weights[idx-1] = 0 + end + weights = [max(0, 1 - sum(weights)), weights...] + # TODO: objects that are already involved in some events should not be involved in other event types as well + weights_normalized = weights ./ sum(weights) + println(weights_normalized) + + # Draw born event + events = vcat(NoEvent, [repeat([event_type], length(object_pairs)) for event_type in event_concepts]...) + + start_event_idx = @trace(categorical(weights_normalized), :start_event_idx) + if start_event_idx > 1 + push!(state.active_events, start_event_idx) + print("Event started: ") + print(events[start_event_idx]) + println(pair_idx[start_event_idx]) + end - return active_events + # Evaluate clauses for the born event + + updated_latents = @trace(Gen.Switch(map(clause, events)...)(start_event_idx, pair_idx[start_event_idx], state.bullet_state.latents), :event) + bullet_state = setproperties(state.bullet_state; latents = updated_latents) + + # death of old active events + #println([1.0-predicates[idx] for idx in 1:length(predicates)]) + # TODO: think about another death probability function, maybe including event age, or add death predicate + weights = [(idx+1 in state.active_events && idx+1 != start_event_idx) ? (1.0 - predicates[idx]) : 0.0 for idx in 1:length(predicates)] # dying has the opposite chance of being born + #print("Unnormalized death weights: ") + #println(weights) + weights = [max(0, 1-sum(weights)), weights...] # no event at index 1 + weights = weights ./ sum(weights) # normalize + #print("Normalized death weights: ") + #println(weights) + end_event_idx = @trace(categorical(weights), :end_event_idx) + #println(end_event_idx) + if end_event_idx > 1 # nothing happens when no event dies + delete!(state.active_events, end_event_idx) + print("Ended event: ") + print(events[end_event_idx]) + println(pair_idx[end_event_idx]) + end + + return state.active_events, bullet_state end """ for one object, observe the noisy position in every dimension """ -@gen (static) function observe_position(k::RigidBodyState, noise::Float64) +@gen function observe_position(k::RigidBodyState, noise::Float64) pos = k.position # add noise to position obs = @trace(broadcasted_normal(pos, noise), :positions) @@ -215,12 +263,12 @@ end """ for one time step, run event and physics kernel """ -@gen (static) function kernel(t::Int, prev_state::CPState, params::CPParams) +@gen function kernel(t::Int, prev_state::CPState, params::CPParams) # event kernel - active_events = @trace(event_kernel(prev_state, params.event_concepts), :events) + active_events, bullet_state = @trace(event_kernel(prev_state, params.event_concepts), :events) # simulate physics for the next step based on the current state and observe positions - bullet_state::BulletState = PhySMC.step(params.sim, prev_state.bullet_state) + bullet_state::BulletState = PhySMC.step(params.sim, bullet_state) obs = @trace(Gen.Map(observe_position)(bullet_state.kinematics, Fill(params.obs_noise, params.n_objects)), :observe) next_state = CPState(bullet_state, active_events) @@ -230,7 +278,7 @@ end """ generate physical scene with changepoints in the belief state """ -@gen (static) function cp_model(t::Int, params::CPParams) +@gen function cp_model(t::Int, params::CPParams) # initalize the kinematic and event state init_state = @trace(cp_prior(params), :prior) diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 53080d7..5fb2b88 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -1,9 +1,10 @@ +using Revise using Gen using GalileoEvents mass_ratio = 2.0 obj_frictions = (0.3, 0.3) -obj_positions = (0.5, 1.5) +obj_positions = (0.5, 1.2) mprior = MaterialPrior([unknown_material]) pprior = PhysPrior((3.0, 10.0), # mass @@ -16,16 +17,35 @@ t = 120 function forward_test() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) event_concepts = Type{<:EventRelation}[Collision] - cp_params = CPParams(client, [a,b, b], mprior, pprior, event_concepts, obs_noise) - trace, _ = Gen.generate(cp_model, (t, cp_params)) + cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) + addr = :prior => :objects => 1 => :mass + cm = Gen.choicemap(addr => 30) + trace, _ = Gen.generate(cp_model, (t, cp_params), cm) + #display(get_choices(trace)) +end + +# constrained generation +function constrained_test() + client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) + event_concepts = Type{<:EventRelation}[Collision] + cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) + + addr = 10 => :events => :start_event_idx + cm = Gen.choicemap(addr => 1) + trace, _ = Gen.generate(cp_model, (t, cp_params), cm) display(get_choices(trace)) end +# gen regenerate + function update_test() + t = 120 + client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) - mc_params = MCParams(client, [a,b], mprior, pprior, obs_noise) - trace, _ = Gen.generate(mc_gm, (t, mc_params)) + event_concepts = Type{<:EventRelation}[Collision] + cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) + trace, _ = Gen.generate(cp_model, (t, cp_params)) addr = :prior => :objects => 1 => :mass cm = Gen.choicemap(addr => trace[addr] + 3) @@ -41,5 +61,68 @@ function update_test() end -forward_test() +function update_test_2() + t = 120 + + client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) + event_concepts = Type{<:EventRelation}[Collision] + cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) + trace, _ = Gen.generate(cp_model, (t, cp_params)) + + addr = :prior => :objects => 1 => :mass + cm = Gen.choicemap(addr => trace[addr] + 3) + trace2, _ = Gen.update(trace, cm) + + # compare final positions + t=120 + pos1 = Vector(get_retval(trace)[t].bullet_state.kinematics[1].position) + pos2 = Vector(get_retval(trace2)[t].bullet_state.kinematics[1].position) + @assert pos1 != pos2 + + # TODO: more tests that include events + + return trace, trace2 +end + + +# test switch combinator in terms of gen's reaction to proposed changes + +# toy model for dealing with complexing +# random walk with 2 delta functions (gaussian vs uniform) chosen by switch +# initial trace is changed by a mh proposal for switch index +# static first, unfold complexity second step + +@gen function function1() + v ~ normal(0., 1.) +end + +@gen function function2() + v ~ uniform(-1., 1.) +end + +@gen function choose_delta(i::Int64) + function_idx = @trace(categorical([0.5, 0.5]), (:function, i)) + delta ~ Gen.Switch(function1, function2)(function_idx) +end + +@gen function switch_model() + n = @trace(poisson(5), :n) + total = 0. + for i=1:n + function_idx = @trace(categorical([0.5, 0.5]), (:function, i)) + total += @trace(Gen.Switch(function1, function2)(function_idx), (:x, i)) + end + @trace(normal(total, 1.), :y) + total +end + +function switch_test() + trace, _ = Gen.generate(switch_model, ()) + display(get_choices(trace)) +end + +switch_test() + +#forward_test() +#constrained_test() #update_test() \ No newline at end of file From 9d88fa9204125257e75394bda3eaab28d44d3c61 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Fri, 20 Oct 2023 15:02:42 -0400 Subject: [PATCH 08/18] write simple switch test --- test/gms/cp_gm.jl | 45 ++++++++++++++++++++++++++++----------------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 5fb2b88..d44dd20 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -1,4 +1,3 @@ -using Revise using Gen using GalileoEvents @@ -100,29 +99,41 @@ end v ~ uniform(-1., 1.) end -@gen function choose_delta(i::Int64) - function_idx = @trace(categorical([0.5, 0.5]), (:function, i)) - delta ~ Gen.Switch(function1, function2)(function_idx) -end +@gen function switch_model_static() + function_idx = @trace(categorical([0.5, 0.5]), :function) -@gen function switch_model() - n = @trace(poisson(5), :n) - total = 0. - for i=1:n - function_idx = @trace(categorical([0.5, 0.5]), (:function, i)) - total += @trace(Gen.Switch(function1, function2)(function_idx), (:x, i)) - end - @trace(normal(total, 1.), :y) - total + x = @trace(Gen.Switch(function1, function2)(function_idx), :x) + + y = @trace(normal(x, 1.), :y) end -function switch_test() - trace, _ = Gen.generate(switch_model, ()) +function switch_test_static() + # unconstrained generation + trace, _ = Gen.generate(switch_model_static, ()) display(get_choices(trace)) + + # constrained generation + cm = Gen.choicemap(:function => 1) + trace2, _ = Gen.generate(switch_model_static, (), cm) + display(get_choices(trace2)) + + # update trace + trace3, _ = Gen.update(trace, cm) + display(get_choices(trace3)) +end + + + +@gen function switch_model_unfold() + function_idx = @trace(categorical([0.5, 0.5]), :function) + + x = @trace(Gen.Switch(function1, function2)(function_idx), :x) + + y = @trace(normal(x, 1.), :y) end -switch_test() +switch_test_static() #forward_test() #constrained_test() #update_test() \ No newline at end of file From 49caf3d40902aa09ba4d87209c6df754b109f853 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Sun, 22 Oct 2023 20:40:34 -0400 Subject: [PATCH 09/18] update distance function --- src/gms/cp_gm_pb.jl | 8 ++++++-- test/gms/cp_gm.jl | 9 +++++---- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index f673c07..f38e81c 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -1,3 +1,4 @@ +using Revise export CPParams, CPState, @@ -112,9 +113,12 @@ end Bernoulli weight that event relation holds """ function predicate(t::Type{Collision}, a::RigidBodyState, b::RigidBodyState) + println(a) d = norm(Vector(a.position)-Vector(b.position)) # l2 distance - @show d - clamp(1.0 - exp(-d), 0., 1.) + + #@show d + + clamp(exp(-20d), 0., 1.) end # TODO: surface distances, use pybullet (contact maybe not helpful) diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index d44dd20..d5af53f 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -1,3 +1,4 @@ +using Revise using Gen using GalileoEvents @@ -19,7 +20,8 @@ function forward_test() cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) addr = :prior => :objects => 1 => :mass cm = Gen.choicemap(addr => 30) - trace, _ = Gen.generate(cp_model, (t, cp_params), cm) + trace, _ = Gen.generate(cp_model, (t, cp_params), cm); + println("") #display(get_choices(trace)) end @@ -123,7 +125,6 @@ function switch_test_static() end - @gen function switch_model_unfold() function_idx = @trace(categorical([0.5, 0.5]), :function) @@ -133,7 +134,7 @@ end end -switch_test_static() -#forward_test() +#switch_test_static() +forward_test() #constrained_test() #update_test() \ No newline at end of file From 0667f3c2cfcc1301bab52e8ee6967e0e7fd1e6d6 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Tue, 24 Oct 2023 17:34:09 -0400 Subject: [PATCH 10/18] add event visualization --- Manifest.toml | 2 +- Project.toml | 1 + src/gms/cp_gm_pb.jl | 37 ++++++++++++++++++------------------- test/gms/cp_gm.jl | 36 +++++++++++++++++++++++++++++++++--- 4 files changed, 53 insertions(+), 23 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 7bcc38b..fbd4586 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.9.3" manifest_format = "2.0" -project_hash = "6c92244b6b2c576cb0a17ab8a9dd1f5a2f86d403" +project_hash = "87c0603f65be5921938c725fd12195c48bc0f93f" [[deps.Accessors]] deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Test"] diff --git a/Project.toml b/Project.toml index 0ce58f8..95a084f 100644 --- a/Project.toml +++ b/Project.toml @@ -15,5 +15,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" PhyBullet = "63daae69-5b14-439d-ac6f-096429ca839b" PhySMC = "79c1e2f5-7911-41a0-b248-4858717ddd79" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index f38e81c..47e3451 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -103,7 +103,7 @@ initializes belief about all objects and events bullet_state = setproperties(params.template; latents = new_latents) # initialize the event state - active_events = Set{EventRelation}() + active_events = Set{Int64}() init_state = CPState(bullet_state, active_events) return init_state @@ -113,12 +113,11 @@ end Bernoulli weight that event relation holds """ function predicate(t::Type{Collision}, a::RigidBodyState, b::RigidBodyState) - println(a) + #println(a) d = norm(Vector(a.position)-Vector(b.position)) # l2 distance #@show d - - clamp(exp(-20d), 0., 1.) + clamp(exp(-10d), 0., 1.) end # TODO: surface distances, use pybullet (contact maybe not helpful) @@ -194,7 +193,7 @@ end iterate over event concepts and evaluate predicates to active/deactive """ @gen function event_kernel(state::CPState, event_concepts::Vector{Type{<:EventRelation}}) - println("-------") + #println("-------") object_pairs = collect(combinations(state.bullet_state.kinematics, 2)) pair_idx = repeat(collect(combinations(1:length(state.bullet_state.kinematics), 2)), length(event_concepts)) pair_idx = [[0,0], pair_idx...] # for no event @@ -202,29 +201,29 @@ iterate over event concepts and evaluate predicates to active/deactive # map possible events to weight vector for birth decision using the predicates predicates = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs] - print("Predicates: ") - println(predicates) + #@show predicates # randomly draw a starting and an ending event # birth weights = copy(predicates) - for idx in state.active_events # active events should not be born again + active_events = copy(state.active_events) + for idx in active_events # active events should not be born again weights[idx-1] = 0 end weights = [max(0, 1 - sum(weights)), weights...] # TODO: objects that are already involved in some events should not be involved in other event types as well weights_normalized = weights ./ sum(weights) - println(weights_normalized) + # println(weights_normalized) # Draw born event events = vcat(NoEvent, [repeat([event_type], length(object_pairs)) for event_type in event_concepts]...) start_event_idx = @trace(categorical(weights_normalized), :start_event_idx) if start_event_idx > 1 - push!(state.active_events, start_event_idx) - print("Event started: ") - print(events[start_event_idx]) - println(pair_idx[start_event_idx]) + push!(active_events, start_event_idx) + #print("Event started: ") + #print(events[start_event_idx]) + #println(pair_idx[start_event_idx]) end # Evaluate clauses for the born event @@ -235,7 +234,7 @@ iterate over event concepts and evaluate predicates to active/deactive # death of old active events #println([1.0-predicates[idx] for idx in 1:length(predicates)]) # TODO: think about another death probability function, maybe including event age, or add death predicate - weights = [(idx+1 in state.active_events && idx+1 != start_event_idx) ? (1.0 - predicates[idx]) : 0.0 for idx in 1:length(predicates)] # dying has the opposite chance of being born + weights = [(idx+1 in active_events && idx+1 != start_event_idx) ? (1.0 - predicates[idx]) : 0.0 for idx in 1:length(predicates)] # dying has the opposite chance of being born #print("Unnormalized death weights: ") #println(weights) weights = [max(0, 1-sum(weights)), weights...] # no event at index 1 @@ -245,13 +244,13 @@ iterate over event concepts and evaluate predicates to active/deactive end_event_idx = @trace(categorical(weights), :end_event_idx) #println(end_event_idx) if end_event_idx > 1 # nothing happens when no event dies - delete!(state.active_events, end_event_idx) - print("Ended event: ") - print(events[end_event_idx]) - println(pair_idx[end_event_idx]) + delete!(active_events, end_event_idx) + #print("Ended event: ") + #print(events[end_event_idx]) + #println(pair_idx[end_event_idx]) end - return state.active_events, bullet_state + return active_events, bullet_state end """ diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index d5af53f..4005a6a 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -1,6 +1,7 @@ using Revise using Gen using GalileoEvents +using Plots mass_ratio = 2.0 obj_frictions = (0.3, 0.3) @@ -25,6 +26,33 @@ function forward_test() #display(get_choices(trace)) end +function visualize_active_events() + client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) + event_concepts = Type{<:EventRelation}[Collision] + cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) + + cm = Gen.choicemap() + cm[:prior => :objects => 1 => :mass] = 10 + cm[:prior => :objects => 2 => :mass] = 10 + + cumulative_active_events = zeros(Int, t) + cumulative_started_events = zeros(Int, t) + cumulative_ended_events = zeros(Int, t) + num_traces = 1000 + for i in 1:num_traces + trace, _ = Gen.generate(cp_model, (t, cp_params), cm); + cumulative_active_events .+= [length(trace[:kernel=>j=>:events][1]) for j in 1:t] + cumulative_started_events .+= [Int(trace[:kernel=>j=>:events=>:start_event_idx]>1) for j in 1:t] + cumulative_ended_events .+= [Int(trace[:kernel=>j=>:events=>:end_event_idx]>1) for j in 1:t] + end + + for (var, filename) in zip([cumulative_active_events, cumulative_started_events, cumulative_ended_events], + ["active_events.png", "started_events.png", "ended_events.png"]) + plt = plot(1:t, var ./ num_traces, xlabel="t", ylabel="average events", legend=false, yrotation=90) + savefig(plt, filename) + end +end + # constrained generation function constrained_test() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) @@ -134,7 +162,9 @@ end end -#switch_test_static() -forward_test() + +#forward_test() +visualize_active_events() #constrained_test() -#update_test() \ No newline at end of file +#update_test() +#switch_test_static() \ No newline at end of file From e3bfd5b342b40b0857c6c94f2d85080c1da8ba9f Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 1 Nov 2023 11:55:27 -0400 Subject: [PATCH 11/18] improve event plotting and detection --- src/gms/cp_gm_pb.jl | 62 ++++++++++----------------------- test/gms/cp_gm.jl | 55 +++++++++++++++++++++-------- test/gms/plots/events.png | Bin 0 -> 13052 bytes test/gms/plots/x_positions.png | Bin 0 -> 15459 bytes 4 files changed, 59 insertions(+), 58 deletions(-) create mode 100644 test/gms/plots/events.png create mode 100644 test/gms/plots/x_positions.png diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 47e3451..defe8fc 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -38,12 +38,14 @@ struct CPParams <: GMParams template::BulletState n_objects::Int64 obs_noise::Float64 + death_factor::Float64 end function CPParams(client::Int64, objs::Vector{Int64}, mprior::MaterialPrior, pprior::PhysPrior, event_concepts::Vector{Type{<:EventRelation}}, - obs_noise::Float64=0.) + obs_noise::Float64=0., + death_factor=10.) # configure simulator with the provided # client id sim = BulletSim(;client=client) @@ -54,7 +56,7 @@ function CPParams(client::Int64, objs::Vector{Int64}, # Note: alternative latents will be suggested by the `prior` template = BulletState(sim, rigid_bodies) - CPParams(mprior, pprior, event_concepts, sim, template, length(objs), obs_noise) + CPParams(mprior, pprior, event_concepts, sim, template, length(objs), obs_noise, death_factor) end """ @@ -113,11 +115,11 @@ end Bernoulli weight that event relation holds """ function predicate(t::Type{Collision}, a::RigidBodyState, b::RigidBodyState) - #println(a) - d = norm(Vector(a.position)-Vector(b.position)) # l2 distance - - #@show d - clamp(exp(-10d), 0., 1.) + if norm(Vector(a.linear_vel)-Vector(b.linear_vel)) < 0.01 + return 0 + end + d = norm(Vector(a.position)-Vector(b.position))-0.175 # l2 distance + clamp(exp(-15d), 0., 1.) end # TODO: surface distances, use pybullet (contact maybe not helpful) @@ -173,9 +175,6 @@ function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelati end end -## TODO: add Pkg revise - - # map mcmc kernel (link: https://www.gen.dev/docs/stable/ref/mcmc/#Composing-Stationary-Kernels) # set of proposal functions @@ -184,27 +183,18 @@ end # e.g. another event type for the same pair or another event for one event type -# cognition: computation around events -# modelling: changepoint model with Switch combinator -# inference: change types of events to help inference being faster - - """ iterate over event concepts and evaluate predicates to active/deactive """ -@gen function event_kernel(state::CPState, event_concepts::Vector{Type{<:EventRelation}}) - #println("-------") +@gen function event_kernel(state::CPState, params::CPParams) object_pairs = collect(combinations(state.bullet_state.kinematics, 2)) - pair_idx = repeat(collect(combinations(1:length(state.bullet_state.kinematics), 2)), length(event_concepts)) + pair_idx = repeat(collect(combinations(1:length(state.bullet_state.kinematics), 2)), length(params.event_concepts)) pair_idx = [[0,0], pair_idx...] # for no event - #println(pair_idx) # map possible events to weight vector for birth decision using the predicates - predicates = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs] - #@show predicates + predicates = [predicate(event_type, a, b) for event_type in params.event_concepts for (a, b) in object_pairs] - # randomly draw a starting and an ending event - # birth + # birth of one or no random event weights = copy(predicates) active_events = copy(state.active_events) for idx in active_events # active events should not be born again @@ -213,41 +203,27 @@ iterate over event concepts and evaluate predicates to active/deactive weights = [max(0, 1 - sum(weights)), weights...] # TODO: objects that are already involved in some events should not be involved in other event types as well weights_normalized = weights ./ sum(weights) - # println(weights_normalized) # Draw born event - events = vcat(NoEvent, [repeat([event_type], length(object_pairs)) for event_type in event_concepts]...) + events = vcat(NoEvent, [repeat([event_type], length(object_pairs)) for event_type in params.event_concepts]...) start_event_idx = @trace(categorical(weights_normalized), :start_event_idx) if start_event_idx > 1 push!(active_events, start_event_idx) - #print("Event started: ") - #print(events[start_event_idx]) - #println(pair_idx[start_event_idx]) end - # Evaluate clauses for the born event - updated_latents = @trace(Gen.Switch(map(clause, events)...)(start_event_idx, pair_idx[start_event_idx], state.bullet_state.latents), :event) bullet_state = setproperties(state.bullet_state; latents = updated_latents) - # death of old active events - #println([1.0-predicates[idx] for idx in 1:length(predicates)]) - # TODO: think about another death probability function, maybe including event age, or add death predicate - weights = [(idx+1 in active_events && idx+1 != start_event_idx) ? (1.0 - predicates[idx]) : 0.0 for idx in 1:length(predicates)] # dying has the opposite chance of being born - #print("Unnormalized death weights: ") - #println(weights) + # death of one or no active event + weights = [(idx+1 in active_events && idx+1 != start_event_idx) ? max(1. - predicates[idx] * params.death_factor, 0.) : 0.0 for idx in 1:length(predicates)] # dying has a much lower chance of being born + weights = [max(0, 1-sum(weights)), weights...] # no event at index 1 weights = weights ./ sum(weights) # normalize - #print("Normalized death weights: ") - #println(weights) + end_event_idx = @trace(categorical(weights), :end_event_idx) - #println(end_event_idx) if end_event_idx > 1 # nothing happens when no event dies delete!(active_events, end_event_idx) - #print("Ended event: ") - #print(events[end_event_idx]) - #println(pair_idx[end_event_idx]) end return active_events, bullet_state @@ -268,7 +244,7 @@ for one time step, run event and physics kernel """ @gen function kernel(t::Int, prev_state::CPState, params::CPParams) # event kernel - active_events, bullet_state = @trace(event_kernel(prev_state, params.event_concepts), :events) + active_events, bullet_state = @trace(event_kernel(prev_state, params), :events) # simulate physics for the next step based on the current state and observe positions bullet_state::BulletState = PhySMC.step(params.sim, bullet_state) diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 4005a6a..8d84be8 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -2,6 +2,7 @@ using Revise using Gen using GalileoEvents using Plots +ENV["GKSwstype"]="160" # fixes some plotting warnings mass_ratio = 2.0 obj_frictions = (0.3, 0.3) @@ -26,31 +27,55 @@ function forward_test() #display(get_choices(trace)) end +function add_rectangle!(plt, xstart, xend, y; height=0.8, color=:blue) + xvals = [xstart, xend, xend, xstart, xstart] + yvals = [y, y, y+height, y+height, y] + plot!(plt, xvals, yvals, fill=true, seriestype=:shape, fillcolor=color, linecolor=color) +end + +get_x2(trace, t) = get_retval(trace)[t].bullet_state.kinematics[2].position[1] + function visualize_active_events() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) event_concepts = Type{<:EventRelation}[Collision] cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) cm = Gen.choicemap() - cm[:prior => :objects => 1 => :mass] = 10 - cm[:prior => :objects => 2 => :mass] = 10 - - cumulative_active_events = zeros(Int, t) - cumulative_started_events = zeros(Int, t) - cumulative_ended_events = zeros(Int, t) - num_traces = 1000 + cm[:prior => :objects => 1 => :mass] = 2 + cm[:prior => :objects => 2 => :mass] = 1 + cm[:prior => :objects => 1 => :friction] = 0.5 + cm[:prior => :objects => 2 => :friction] = 1.2 + cm[:prior => :objects => 1 => :restitution] = 0.2 + cm[:prior => :objects => 2 => :restitution] = 0.2 + + num_traces = 50 + plt = plot(legend=false, xlim=(0, t), ylim=(1, num_traces+1), yrotation=90, ylabel="Trace", yticks=false, xlabel="Time step") + collision_t = nothing for i in 1:num_traces + if i % 10 == 0 + @show i + end trace, _ = Gen.generate(cp_model, (t, cp_params), cm); - cumulative_active_events .+= [length(trace[:kernel=>j=>:events][1]) for j in 1:t] - cumulative_started_events .+= [Int(trace[:kernel=>j=>:events=>:start_event_idx]>1) for j in 1:t] - cumulative_ended_events .+= [Int(trace[:kernel=>j=>:events=>:end_event_idx]>1) for j in 1:t] - end - for (var, filename) in zip([cumulative_active_events, cumulative_started_events, cumulative_ended_events], - ["active_events.png", "started_events.png", "ended_events.png"]) - plt = plot(1:t, var ./ num_traces, xlabel="t", ylabel="average events", legend=false, yrotation=90) - savefig(plt, filename) + start = nothing + first_x = i==1 ? get_x2(trace, 1) : nothing # only look for collision in first trace + for j in 1:t + if trace[:kernel=>j=>:events=>:start_event_idx]==2 + start = j + end + if trace[:kernel=>j=>:events=>:end_event_idx]==2 + finish = j + add_rectangle!(plt, start, finish, i) + end + if first_x !== nothing && abs(first_x - get_x2(trace, j)) > 0.001 + collision_t = j + first_x = nothing + end + end + end + vline!(plt, [collision_t], linecolor=:red, linewidth=2, label="Vertical Line") + savefig(plt, "test/gms/plots/events.png") end # constrained generation diff --git a/test/gms/plots/events.png b/test/gms/plots/events.png new file mode 100644 index 0000000000000000000000000000000000000000..0e70bbb710dbc8ed9c5272101797a2bac552c579 GIT binary patch literal 13052 zcmai5bzD?yw;npB8>OUE=|*s90qJh3p-VahgrTGbNhPENM7mo91VtE-5D-B+rSrby zIrsa{Ip24Gcl(FG(LH;=?|RpI*0Y|qCiUu9wV1ct*s-^Q}~+v{Y6&%Mt+q*HPMNQle9{UcthcEB^7)* zCMGRGOh#g24w~>pIGikARtXA?dC$Pi#3XM`L_(aDe-}F%t~5#&0fiRY|8e^y*7sPx zb^}Ahwf(U(aZNI3;jD38-`lc8%mF*9>4{mfzgAedUb+is#rM#Qvn&(FQU0h(e=oQq zuGMs5wK3cDOW!;Q(d212nL9CQ$W0+jgpV0M{?s=oCMITK!Qwb)-oCfq+3n0YzkelG zd*X7fcwWlKcbaRP>!mmPE87o5-VfsR*5z8CwS-s%>e$%WJbai4-hV4=m=!xry|`u1 zvh}0V&v!L&?0ds} z_9p+-wQIL4>!*ov{Hv>xjy{s4q3Rm^f7(S@E?KOC!d9TW8+;=JLW++W;b(3)w+^%S z^Cs~_=Mzc}-)Zm3l^XGxa`HW~Txt0d3~9fwf@@3QW$sY%sLm)hW?c2Ts9t#ZT~|6C zPGv`?WDZQmfiCmwm`$e=W8Y#KlTLdEK0+>nrs6lkO+Tl}&>!T+m9Uq_FeVFt8=tUu zTRH5NzkxtxL#H|!QM33wR0(^ca=9P^>e9u)Ol{Ock^l9;vH=yMo~JP3 zAVN7#Gbw1i>c!5tvPlYV5td`3e>tjTijoxg4^O85oXx7i1iSV4IV1J@Bw=5l8a>4x z?ZcS2cf6wALezMeWAarnI!7G^*|DEINv41@o8!wRzuM`jG9DS~sdXSiXN#K3Q+rry zBpH%|^|b@lEXW^KS9cXLbehCTezix4Prv=@4}2}; z$|DXKA&@q@{%2reMxV~FNLW3D{A2+4@hB|% zaQi)8VL)PH-A_bta2-|fwkP7KG#hT-jQi0~d~@?G489a2t2F4b&sl8 zbf@MfPb|eTU~{Ia1bJTw*3nG=G<{=f3EyoW5QMijU-Tw0TDcdx8U{K(3Q$XI%Lu-vJ_s?dcG0{HSVlpNfgGK!lwrCkS zo*8XX1o!1`kS1MvkGEvQ&br>oE`~MHk*HKZC*KZ2qh|*EX9hDdTl%Fh4e><8#7(4sMN(loE%7$Gu3#-n_gk?~wKqxpfS0E(ZyMf| z**4e?Br$ywe7W0eq3r1K-l`fZU+S~z*4|*P=(V#w1qIzwIR_BlH1tBkMEU3Uz zC-;kZH^jqUZjBEJ!&1k+i8vr&^YUbt?PqTzE~(LDF?8pS1;aXIpqc=eQlI45fTZ-d zTvJb+{0<$nd!F8i2Xcad*z)BBbvs3NF?bV02a)(K5klwMj*~T9O+yaeYp2lVozz%#_>hmz-S3Tn8BnM}L7iJI^N;*&q# zpffq%L=P2meCv{9-E^+JK-K)XPwE>B$CiSVuxyc5`MO%4pr6OLIpW*9CqGDZct$C2WV-Y0TWlCZ6r z9KxqYrp{m1yo`9N9x^rjuu)iY+h@kbwGi~*zU3`W74Jbsn*E`;IzF-FSyV@8soy2B zABq|LuNg=4mc~&kN+mUV8oaO@BqM@v&5SI41%^~KG51^pJ;Id?sUn!MK1qfj&pGau z=azJPb~;vB?sb@&YYDJ0i_Y`LM4oADpkwC0grK0%mBlD&BJYnPENE=(t9mEHuBrQU zkyD>?33XCZ&D`u+8RwsY{{gFn=GLP~_%wj6zp|E7h-_jrD$}Wz9d&XzwEz*t#G?Ko zCxn1WjC#z!>DB9cVszOGy;7ZLGP>bmL^Pbv^3ZpE49sE~(znKnCo3W&^IKcj^t~}( zLohH5^cUK{7^#g}T+Scn?omnKoM{NCrD3A}QK6aq3Y#s=s6o!RlawnbYnKhi8tMG`H^D%PfOGik4O4KUhkM z-_8`iVIGK*Ljakcw8*~PizcPxsyMA|5-u+T(5Q$LJsXA$dBN^SVFTz8TgW}N8=Nra zic)HPOapE$N?B%Yr`vyn>HiIRt*%dmAf8rv?%nVbM_(`EE$6TmSoz8F2lA=%S&p;! zdra#s;1w2G9yucc90AiUN5JZSbaf3AZhX?6aEOCjfx%rt8r#;Vs#TwwsDJv*;eQ_(oOzFDF4M-)-EoHe)4$} zPZG;l2iv`f!7Vw}=Sm8=1Jc;wU}D;GpngJ+63W?8-!;C=M$5tW)*^BB@K2fwz4B4} z&2oU=dB4$nU{>G!d;9rKjn5*xIybaJFVCwbJJeh)>lU{f?pn%aseifWU8hz|n)6%@ zTG_V|_eh=W)nq~I>;Uh1HmQ(MiPOv`>60f!eZo0U;LN?YSbt5#$YkRj!R+ZTRaQ<+ zw&bLL;~`obovd?KQ2uL7|Ej!>FgQ&Q047>762%1ZdGaF-3y!mcu8X+(_94Cwi-F|& z+(PMC`SQ_{tkcbf2%Oc4;j|s~LUzC{#Oi2Fn+4zTp+b(fdghM<-oE_#Md+7kJm!-p zz6Y8O6!)3!KsbwQMhvIZ2QQc!<_iM{|L59dm6m26utr2EnW880S=&+nhL*@1&z7sB z?&2sZcg5)-Q;jWYZXij}nSCK|gwhrABpoL#ekk;NC|Z!Xo!pU@iHZ(&*f@9{6KCc7lP z$$g*S>wK%xt8}f0sS3h^O{mzQbUZmOm(O0N1JOA`z_HOWl#Pgx%uqS9{i@kTj1B9@ z`!1NbM+OKeAt81?*RdIgx|7}e_;i~Fx9&iSH5z_O{^ z#XZ`dG}qL85`)Xu5_fb>tNF-j%orVg;p2;D#6@7naD|wLx_Z_iL;i#%E#p9*%JR$# z8&iYTUka>o<9a8*<8MA&T;$~By1KewD4H&j|i42X$7tkfWY0Ux!3{1?COd-5Z zRlc@^VUhXGnFfvZGPnt=81;)6BQ@)ehJ1VTDDAin$P(kzD6tzsMVo z+SJf-@#*0#5sh5Qsdj$PpC6Zp(kaQw`+9qoj46`JO`0RVwOH1OO$l}4LOL$Oq@3xz z&g*BU{V^c{M+;9WZw_&$UAOGuYQJt=sD*l2PoJEeOs5&UG+}g>Y8L~*re$GySZ?Wc z6dHj6v8cmoS;B$owpGs4#48kx>gegI7!V~5oi6^k(QY7+l+gTo za+Of#!lKTIl$4Z(GZGH}Ft*-KoGTsVe|2#hbat>dQISSC@J4Y0)?fagza{VC1#gE%c`h`9{;vZb`rdY0^{nH2$3I87)wZSFrZ?Z`MH&HJRv=qoU<#)IRDsXUD!1 zF#q|^)m2I5PZ99g#@Cbfva(%?jQ1ZL8BW6if-Gb*+>{RV1zV7neFDA>cJ@_Bc6Dc` zYHUqKRkf%LX@j0hhugCqAsu~Imi?N^q-nT9LU*XLH@_I6aMc=7WqTJTJs*DGWdoNi z+@|dO-Vy`gEm^=o|)o-?%^6S&!%eQy5V+LJ5!0H}v zmJAL%Pjo3O?;iVj%0l;xw!T(dTUl9YYL4#x=)ZmYcFntYIv4j-A2IWj5XzV|h&zzf zJ3Rp$6jLA)rK!nxh&HQF==w~*rV?61J(O;}ZQ?U|pV@eNe0&^_l9!p8`9Y4j_gt&r z)^uG8FuVwVmZrpM?b`Zt{u(a~p|Yls1RDICk&z{8wY0)ciOrZyOiTkMS<(@J{1EfI5#8fdQPikpB@D3h|=Xk@EsK=%J55Beq92cpkLmt`LWexw#$e=1M2bav-L1ympjPt<2fR;Yk#}Xr~Lhl@*Z$A7&fQ^yy`;=YYw`Omy zD1+NfOv|59S3ke%nZG&`Ddl2Jj`We1w&a)H2=cTd8Xy5WMbn~H>v8U{XKZN43}_%P ze0)O3SI!Y;!NaM=XSf*X>MtHiaOebF0}!JZ=Vy-GMoW^|*X#^vABK8j?t^Z_kNuGu& zRzCdgozhsqwerp|?u~m$`7eBr%Gptofo$I2nae=T*9ZxP|E4PO!_lu1R$Ztoe1cTF z=^dg!UohqW{Cu=IRnzqvSlYwNEIFP&xN{wG=SB7=F4lWgEP@4R*FdsYX(5H{UcMLZ zz3-drshGsZXz*`fAr^=nnFajYeleu71XD7@IGokhpN$KFs8xRNhlkYDEr z&A+OK@^+aOFnf-UzMWd|f@9~9zr9!fv?FdBadAmZgsfYnp8|d-Vo$`$wV|h9+1pPq z@mwMYQ~l+Atl3eud{+$(;b6d z-dD6)LCGjZKsc*J3lZ-_MHbO{&w-06nFm?_{`&Oebfv128}kYsVkX*Jsj__XX~Z5{ zWfo6&nh;`f;+Bx6z$^mzNHdqd+hWYPyy)nx+r)!EPn7h_5$LiKDJ{j!UEquNE zMh5ZvS7XF*q{I$X@40ndbF)J6m_uEw-TkGtGerX+!O@4P@SPe zDNw&SNUrDp3tMk5dyP3=e|b^n9C`yL)cQ&_hCQG9wi9X6#$D_IHk(HaMd|B$!?slM zJkkTV=FtK^SJ6}0Hdek9UL8d%t~@dqdZ8npqf+3or|c!}wYo6f41e@YR!O?lRN`VU z!MprNzxB{_iV#`+w-)F?xose5suxecQ@d@0K(KSlJt}nkGqj%!px02_dJdWCqCwso zSR9;9v1RfxKq#m$J>im_c+(0XtrY^^`!9t7)Qd~^3Z8U?#7Yby3-{GS>?2PNLoJXK zt*M|q;8rn%cjy=xJi)pNpPqPA3E}3vMhWc;D^N5wGvWfWaJ+)H#Ce2XES*&5p`46N zZ*ODu9i~%DqGV9tGjfzJ^GOgur2Jp@Gef@C;#9Tcw$0TQ)n!tYoZH$7=9gd-6O&AU zj09W&E9*vZFgsUDGeMN&wkr)jK;5zkrJr97L6ZHu-c8Gx6s~}lg2TwDHTj)k!j<*S z8z>j%0&D(H3!-EVPL8hPy!u0C2eDPWrDaAt3`8o&|dIXmlVgi8z@INImU_00+nbs2&< z0GERb2OdqVd1Wfl(Jft&KzVx^nDYv?gt*ipODopHG#L-nAIqct2)zirkeIJS2&q&l z8qr`9x}1GjhiMj*8jJTZ5!aytz4@xu*0zX3+2qTF~&_7Lz>2g=9-cW8hiTH8-fS!JY?ux3%2YWISL(+p;~>~ubS36$xB?;I%!Zn|;L-+(lF;b8G91I4Y< zq5kuRrr`is4_-Yn$tgsrN{6Z8`+rm@lYdkwgC3KnjlMt@xMFgnb9Zm7YXc<2x@Jm; zQ`fm_G?54$OV#)fB{woT+3nP0Eb3AZ17c$%y77~m%}AC^_b`tMA3uUY$$%g!8(GIZ zqkGCO9F#!}d8qr;%JqeNe!F9C5g(r-;HRkn&J#4WRzNoaNbvW%=_r-agB~e}`)}SQ zqZA1;*g;)TWM;GWcSY7sVGJcYoQH^4h_Jt?Ce%i4G|S6s^v|&$VB+HL=AkyEk9%Zd zZgA<1^AOTTzTr;LyU%RhDx|`%qsi!v;WClE(*hlT@ZSVAgH_g6Hri>E*eqHhB^^a@!QVPF|Dar3M{z;>3GWHr| zuD_ZD#2i4fRQaN?E=UCh-CO1@ zWofDJs&GQugo$3cZiU|Dh%_&LI^%XXb!`WLUQUw&CoLv6ct}5*ly!Xwg zbKm1*Uc#;$&IThg4G$y^Pppb`pOv9NSV=20W=$SC{(llm&03Fm722G9mDx%7v6M6;c9@Z=$T*MS7 zqoehwC}@QoaZ16Uu)800F1*ut2}#UL(3NSb$<%`#fWFkM2cCzglp+3auUiYWI*P{p zo!OI#g(dzAR2BV-o?=@kMMnqJ8$Brw5GKRNJMWL@C7xq#%VK^mCj4oKdHx zbcsFwy#IxHpX@)SZluPMlmz5lNdy8LXRrOWDwCL1tfb#G7Ci%3xr-1}_fEWse=%xn90-5@r-L{H;&RD&^Mu3q9T{GuACX%*F4nbF6kRux=u~ zLoFIe*U5f6Evq$-{vYwc)KE z9arnOPxy{bYWJJZv*l3Em$*RA#ohgSE|RXo^1lDq9uAWR$U? zo8Qv8yOLS-^(oLoVILo?)X;d#VBEZ) z&fD11(mw=?G2)r^FGu!V(15U18Q|kIH~XQzuQShWovc1*-(Gq6m6n6;2-(?jIXYoh4?<=&Qe85UdutLhUgDM8xYV@`xkuu1B2WDe6RSp4- zTCBNEosodRlwsrtQe;?+juQ2~COj_)l|ef_HSqgtn)cP{!$+u|UaYx+sq^TXmrDs8 zPGd!sjQ?VI7+3E6n<8l1=JYhCrnl@ZLB}8eR%hWc1+EYXNfB}u;BOsw07Uzb?g$3O zP2d{5LUb?GZ*1D>r1}vwuZiyz%vJQOps0TsCzn3u9O4otz@1T#0C=2T}B#>TRWJ zE$Xl5<2*b8!C(G(IB3x1VFQD$U!49b;_9c(k@VXKrLPT<*$%Rp%ni9ZZP{*C1tS7O z4fDY*{)g+V;fcePbWFl4O~U#9Cj$4!6^2{ke;sK^8*=Ti7tB03IQabeA#yI^5{M0~ zDWw6&+gg;1rch1b(O7z7(%Zl8DDW#Oa{ca+UqXO)Uaeo1zEv)Bd7$m)<|Zc>R;%krTQKUoyU_In4NJuo zJ~=Qvy!E9>)z`NXOhTkDE?PlOPt`66+=GaCjqABtSmcMJiYs-^%v`Dsm~P)L;_cnY zmJZ4cx?QDjVzRljLr+G=?9DhQPu5ps_lcF275PXppKX(t$6KNg1JG=1#26SXz}tS? zA0)V?`zLMj|D?aa)kp_3|1}oM-0crZx3h$uIr;eP8Y=Hc5!l$-7rQ-2TQeB4@G22lfQ`*^7 zCvIs~vc4bv2`kIXb&iwjRg*I!6XWBc{$JUs>fq(!adz7R>M5Uu1ns?V{p`&nSoDEJ zw{A(EyS5WIzw+u(BtAVb@YJx;>0`38D= z&$F{715b9>#)|uUdv8iSez2}0WjXcE#=lLPiOK5gTg$_Sj{QNdTjE|B5MuEl8 zxU^K1l?}FEp04C+Ja~{TV5bgtn9gHtJHPslq1LDZyyR=EUk!UqUtb?kVignUyb%LI zBR~vBmwn6k3BbZ)Z}idf$B#&%IX2&(|AkBQ+UvuONmadkFHSBl5L;5NbFI@1lD@kZ z_V)J5%F*mG8`p?? z9#6`qUV?y~6Gjd^ih;ww*4QPm$5gc)8QI#t9=0Ww->dGY%Mx;Y-uu+Qw=b5Gi<5KC z_ecDhtbvgc4H419hYuNI7J+l3qM=opG%>|2f)d};W6#sm(++?M zaN^q#sJy(*qeo72ty8*e$6|-B0q^_z5aZ*Jk~eQ=ff5ygLkc-LJ&mUpupN7?rlrS` zWS8q_3;ge}|G`wP!x-2m5c43jo}Z+ECvvj0Csy9Wd1KtLeRaOmV89J*#LmrK>niE7 z^a1P-oWTtSNy)$uiHs`ePjCnbKw*bKARr%YYZ$PP=#$m5gw-GHC0hyxi`Ii~dq-=vHY>I$(t&gvr zZJD+_uHysi1|sCw5F9%@J2p1986V(Jz?-+`+Gu$E=DAFBo98{hwz#-LS7!Rw)4y6e{1mpkqVj;ek5QIr zo$y8x0+x0(kn!MkSy|catO<~2K~@6NtX+2m&ZN2aQ3osy#Bs;bv?Fkr)|K1<@yDAd zKeJrbvqkUn^J@`BOw%*QE8LFBAGS@Y*C`b|VJvNXOcMx=$>Embc0>uOnbTbXT@83^5)eoif0c$~k zr#O8+R6UXgKr}cwC^F}ZcF)Yr%-wx2MMJs?1Q@{C+ltS=CV=j1Nl8iYiUPT#vsK)q*a`Ra1oVV4+y2glo?y;t|66rSBY1+IZl zFFo5{RGlzn)1Uw#;>H4a7Q}-81&E2UCDCGjTj4A=waoEq>wa=jVqRT5jfy(n{wBf2 zg%`>mGc`4(at?O^!3Smxy@3D@QrFf*g|U~{L5|NH4ugPz=VYbnQI*7e&~Z*)Ufv}^ zh6g@1|I*3ckE3Pl{%6$qkznUbOG_Ykf#^gnDJ`9pdY#~S!-)B|YR;;hBZS5L3 z@bTkpdr(V|HGt$V($CMf8-F|lrZZ{w!T~$`4!SX0HPs>~aDD*+7H%yP8DMd#1;j#bg*X8A95w9LiEG6H_ zz`*InuTNX$FJ8R3y1MF#A`IE6X!tpjb3U%Ep=_KRrK_V;SX6X_R+w8%>@;3vE9e=k zH`f-RktNK-$q8HgQUv%U?H!BI@xEqL7nkiAaxPq4TmYZ~uRTE} z>FMr%&Sy?XLYa^2lLw7Z^5xP%&t`*yqpP8g zzJAB?TmYb`8(R(Y(y(zPqk%d&tOe|}>FenU+kbuuyllc{zSZv+pqn$`gUrQuLfng* z()@fY(2oI}B03=V`@k9`(X@7 zbl?5nBNcFzE$nRMJhlMr4xE+?yE~Ob+ofrl!Xf4eF621zwS{8lTgdfgaBy(BUPakF zf9CV&wyqU5Miq{em0)T%PfJ_7_`wUS4(oo04CGu@szfg(ATGaouiO8nzc4Vv%6W|I zH>c}V*E!)OX(u-sHA#`?Jy;!me-sRe0;nF=Kxe5qFL-^RjOHuwQh*PHJL*uL#dop| z(MvES1wuOv-s_`mYHA8O9WsYSCB{cbzYuXXQ&NI^VW4&bQ#%OKP)9&3I$a~6JiDLo z6HZ0B!0x1^Y6C_Hs}mO$ZG5Z>z~0`zBQ~-1gbFC@pyhnD{f!pbmzdZcFu9oo8TZBB63v25I|4#N zbHF6G;vxaNNdPPbWTj9PwC*0{6-98XtTs9%dehf2wWNaAyxm0N=3_Or;WRFNy1?Bo zMfxBB%tZhQO##~mF9&f>!DBQI){7|D@7V3eqyBwCKoUtxN*V?9-+89KufHEm>|Qpl z}+W{K)m+M!lYiKlRER|1gpm3YDA}OfQiw8r`dC-VxXmA$*`3lFp z5fdMO33{C3KHC={fMgEGwQXxo^R6%KoO8h5=EE@X^eT*WA3ogua62K6PD-4Ot$3iS z5bAN8M6*WfWFZ`aE?d%34b}iCcmrR^(UctotXikVrx#8Pc~?`jRSrrLa2lWtJVF5O z)7#&lo}NBQcik^@byrmME3D4!dK+LCy$~#Qta}@5BrY!Q%$S1Lgo=d30~lnQM>U-n zvF4>;X+kW&2bdb>9VC}*JK_O0cHO_vkW78GrnuyPv|F<=*LFUVBjISLfY>+$tn>=( zXsaAFm^nB&0HZCJl4fQu%+Fs0#sv)6!EL&>w$^6o`EZ_0$Y3f5NY%9=x12tLQlrIt ze+iS1kC(Sy8rwW@PhOTt8dwfMh=i!9FW_gY&$*+Y zJ$r0)UsF@l&+qimt{qfHTQDalr>%E@9mT|e(=my?u7gSgMSz*Q&gd2ynPBk62>`1q zv3@Lc^j1)ggNKNLGV^^rZ92E1mVrT|-!I!@q+id1hBz2GY3jn+h{|Wc2Gk7{Ry?Pa zL<1-R@g&M?AQOKa90WAo@k(_Z<_)TCkR04V+6)5tc@!@zD|^VcR-bnYLb=IzZ}EnD z9q_y9WxkG!FILseY-|mMflOjxPDgp|Q^?gX;LU(vwE%}$Utb3k0w5Sqf#Z4m?;@#( zWltnHO`Nb@aE>756c-oI_^`0C(M$Te0sIsb69f5EAQTWL^$I;ePW^%L6bkl$rvm!x zyYK{SJwGX^ zm0%QQt&5V^Btz8Q%ETmFuUrjWN=gQGu(Y5_FgGIXf3SKDoPKSxYW|I34a93}+8e+( zO+2fBfB-u?J1?&bIDZxvmP3t3K#huX{Q*5M3Aw(aiMIu%ZQZ9QL>yR(&lzDTzked3 zP*B=2dUX-t`98DhiGbY*89DjW$VgDDpPilxJ0?^AU!2lhDa58At?a#b`$Aeo2 zdO%qN1_zA1euCNrT^fL`Fms_e-4$Zq~_*cRfSw~v9cGjf9Yv5J3=<>(`X^5Cq#0L9h}~*l_36 zAo2_RhkxgWnlf^R`IFOBkc1%2$aQ6fo8B+hr@XG}8_W=FC52V&KjUZR;^8u+_(=F} zpR$DC*nV`WW4Bt{vBtBov`T5U*rRSmf3+;{rQSZ}%?M%&ovR6@Ycal?$7)B)ADKTS zWK15Po+MAwCuKj*_|vphhN6yOf$#R?ux^C^5OEe177{d$4nvSuN0cIh%!=W{@3&|} z&_Yusv2x6{;hOU?Y=W%x#)N8%vlfyq2!$DTm>G@)}$!rHiV|_Qpa8z`i zQgV_8T>3%-^r@yu+uy!@dw-&~DIZHk_d=xd-CfrRW!5&YH#Eh?#a{&!F6GXo6P)G? zG%!F15Es<{<@^4R?xizM;2~Gt28gfs2y|-&-qh8*MzKVcO+|e>`N@KqPt2r_Ifu3T zN`yKw3wEvJfQsLGonJ{U+fra|r-U*Z`xkNKi~0wbYM7-IHm|3zR*CvYTrjE@udSdm z{+{vfFvHk{j+SN)AfTw_J%|HGY{@d59nqTADZpO)8 z&|jw%oK`jHj^{P*#m^h33>U}@@nW_!k6q;a{V~mPLvutMub1v^9<`ZYHThyT)is(1 ztzp(nj_$)KDfd~y?5{I1jj_veS6-mJAI0+__d}R*2VPA6*fx<5oKgSV>$6AJ>?{D^ zgoMGsh%TzOD{arXXilDl76ft+eJ~)qYBD>48u&r2_*d31zi)K?R;bO+DqXjy6tZ;z zm*r|t@l!oD2`4+VXV#Y;ld~InCl-*aW+=rV%9}gyk5<0li(`G167ER+=rCFe`A(*Q z8qj-M6}263`2AzopY_1L^HE!7I0JeI53#GC^!%n#^=%#sbo@9mH|I2tjg*Hnb`SzWI{*E(Va;X_YH*2;S#%K`&>V$xTcCX?2UksK=+Di3LIf$7z93ZILM zhwn#|F}HPL)jGHx_&q%f&(^pUxG}f%$&()$FodUf59Bzf&)R8w$eKNOms>HnaTRNc zklBz^113-D8bw)8aXs|aB}-M*RzB%u_n!22&GJ{SiMY&fWq3$B@nyvziHD+Rp57P|4iXE+etGeQU^hB98_$qT!SAV9At-XEG+l5!xlj0a;D5q2X6m+Q-p4G6O4&fU*W~4$%z_OtKPmfu5&K+UfbBn%FH}<$S5o<6l^|d z@9!^b^GZ^GeZ04>e`R-ne}8A^=9SlHkJ0U)?>&8H5|%$x5pX3z{q@7dYwQ;)&7O$A zqGA->}uF~U` zV?Gr$Ee*|i#V18YhXfWytATCzM@L8RRK9DgdHlF|!Fm$wFiR(~CN`=RH^NFoQ&Y?$ zxlE}=jO<9$%urI5nwr{st{t!AXEc={K0f}Ky^!yolZ1ps>^R;)#pZ=&W?zN!+o`-b zt^`CGCvJ^6liA$o~~MScUlCV?!qtPdoEq2 zn>R@E*e_gV6Vk<>C8um2f<+n|7dKsbmw;t{zAMS4Y4`b^!^wKLsi~<|(ZBOa+Vu4F zqGqf?JVbYkuDRqENacqjFIX+z59bqKzGP)$%5HCOKcao_9)qSO`Gm5duxKZi<<_PX ztSvkn8r$A}c=_^W63UMAzv>nhT8fA^#0UqtvZr&Jx2c*q@vyF{L9nFt6v^l+V;`RQ zd5qnvaBv4lG|1-jQsqduv!a}z8S)M_;}nDs(NLT`ac<+_E#i8yKLYX2q>bWh(A+$C zaqgxrQr502(7O7jQH|g$VF*J4^e~);Bu}R;5 zpSrGyz3pFRmSk0#<3872x&S7dfaZ}DOsafGo*-ls{un!S$O46y8!Fu zHA=i$0$i5X8_pCWzdMLG%1K(84r+1UTmR)>DM50BpR7wGSM?!437%y3$QH3DP5ki1 zYGpfU$9cbXkFVw(Z-Dc;Ad9V7f`d<8PWc+DoxZa-uv!HV`f*O~{$_A^5Tx@&rVH63 z4nYT@72|A$lSM=t#G20|??cPw5SB2$dL{oqlrUv}Va&M9YK+4V(_`E}jbjDHyAI!< zm5E5`s9WYR6Q;wR<#xixnPKRDrGhQeOEQCpWVOejy<=JCGh}v#oLMypu)K*(+FBhy zRoMjx5kF8N#onR{L3`g@1w5qnPh)bxMFhOG87AJAr<` zFYh2%QSe)?su%SMYT7RG$xH~!iWfCJG~~0tasTe!i{#|b5)!r?LgVHCzPtM=uC}&T z&SyvR!H=*-UR+WWulTI(z}nYJQO)syWDIWd z+MLczOcaLQ^!w+hf`Wo7y-Tjz^)ayO;#PuSv7Cg4&oPUOi(kI%))s=Ns;YYScfW;~ zheyElolYafYr`s=*p`--`uaoILuRI?8`DiChtk1Gp3?=#QbcB&uD}1a!3>*yeIPfq zx%mo9DEzL_<+4_XGd8WFS5JzcKP6k@;)qB4}!($mwsoqv@ zxYRXk{-}6()EM%O4gdUZRsG8}vGUryqp7~`AnT@uk&%(_{V-KUOynsQCMF@FqwhAP zUVc+r+SAjMtPPDhJ31aaRHdh9US=m>66%PN;c%Ht*qqc<9kyc%!-PQiJ7jj zX+9sR{DQdMVrU`ToiaUrrEcK(GO@ld8yW36OhA?Ju(D!GOi3{?Hcr>zWD64_PI~cT zZDHY7LF`Q)2Kj*ZuU{{Wm6}v`;W#`fe|BAf45R=2QDq0o2?2~<>uA^fFYqzZv)z>0$mFY%%?j9h@ol>Co=fpkLN`XrD&)Iae^GEX@ibQ5Ej8PzkLCm-tCcg1m!{TI*bGW}y9ewT++ zZtcM+E+R@8f_`VCBoaQ8js$gn@Wy)a)V%W}h497imL2#(o$YYrEq`p*4-V0fW}`^X zAUdp-;bHp8QBc}SYz-}vjws=~s=1T?v{Losdyr0yh+;u|caQyPhZPVzo*3HJ^LV54 zevAxp!`4Vnq>uxmw&rs!-jgbQuFL@$P(~@vKAB;(Y~ccEhEF;hvWllcq7q^E{$s4gBj;N#>oz4?3V@ePMGr?h`N$#dSI{6f{UW(NBV3Lfo33@*fvMku6oB zkRVc-lP7ysl+|LNuoMxyE6{!gf6-30r{QC)+1Ou|+eQ1EGiM7a7CgMXpW>Nj>Re_F z=x*M=U1m}zASc&!BSqk{P4B?)-2$xvo3sTD#dDX37@%K(gx8Na@<&3TNAb5vE%QEAIzo+ zp~el8=xC4G>Q3hMKibiiB_$$i^x5V4AXGzV+TSWsJy{_kp~B&-oJri->g*Skl$4TP zI4Gk@h>1BhQgcj)D(`;!9t)xL8hi0(wvs<}#NZ-xZ?WZ_=w+N!)(+Rt>==44+koYwAMp_t@#vro1JW}~v9fz5C+#qzAUY zrXjT^ypg4q6?`}zKEim#41jJ$|}{*{Jr zV6J2s9%=X*CB)s9)Wk-!@KlvfytcD4twlx_|6O~c23E6~3<)KkG;>o-q@ayz-B>3%F2WFwh3DQS zShRCIt*G=um=sxoS}~L5b(ag7-3v|}WF=e}&%Luvv(tOnUTVbC}ut$p6WTxFU1bp4l8NajT(CxByF4zgAdjXURVB^c-_aY*pJ=Pdol{E+CW&F;0)t&42gKkTnq(4{y|nk)B=y)z#NG zs%+Nk>$o|1w#4?aPX#OIbrH3VlamuMF>#rBld*06oyer5B+eD=0AEgyPSTcYmp@;-M zFF1NLNJ1w*KR@4_2Cw=;s=G>6#%tqb!llK&R(6&xorNxGdFFT)672x)sDhon{iDAy z*h^n58r+Kt3qN;v7puiH;t-IMM+%86L?82JREjW9t-C^+wNq@O`V1kU0&@vz-p09# zGV+Q~pE#!mK|o1pGRX6~Z$JApp*cpw%W%V-Sfb_cQOPxSkECzk882*UIg^<28*cP+QT%M1?SI?IsLT9TtNxc)u&n)E=mQ+F5z!?T z^wA;a;%5AkDUeOSh+{rk)<}UsUYw?uL?mVRkmh?>bL%mXupQ#E*0V^Z`rRlde)$ zw#oBhNRGWAo)v&a4Gm&UzMF3F=ouK$IeeO%Yvt;?RW}n*vWy+f6rk~H$(=l^!h{?V z-FVF$^Y3C0t9BfBhKMxP=w(~PGn+EW!zulEnmVS06A8jAcBy>+{CWPIbkLWvnYUVH6G^os5?hTJ%(OgUmY(eBPETCiLs_u<v*9`D-w*$!9kzJ{}s#l6~|R z!o^#tc_jS~DqcPXlEKH<7t+OXly3{xrvYMs-A%mYb(DC>*n5;B1WC9$etOOL*DYMx zhJ_4S0$O7YbKjV(2Rsj-CnX~*pHh94rB!mf6s9&yK5+WsE)pFlo6it)EkO>d*Vc4Z zhE#3=0WD3<-QT};d$8Ie!vfk_(3g;?12l3(%OvR<^s+D^A>oMj+Up_z<2~+d&$X|E zpFXJ_KKq+3QI>8YI_4tJKE&if?z%`8mW;`*z>GT*=T^*>B;6vI+gTo-sJ6cx_}5SI z^X}?60SWES$)VfDi%9^4Na-ajgtHkK8UMnlOB1oNVlMUzky+ExUq;8rQ?9Y|%&-dx z2=MSYY)&^VE-v!WeYO^fRNnr?6v&lAOe6B~!w03qn*u34)*`{k#kl;UqM}@`;=Sx5 zj*k>X4|#X-rLQ?U81J}=s$od3cc@og|LvW>|0(`s5eh+la!eLfdQ8fl)t@2pWpMB} zB>S`7aSI3_kIw#{cH*wRkaByaXP+4!_2i+uJ0u}J0X{xHcXxNcU4|b&e?BHVb*pRP zq2lKv{Q|D2mGN>WmXH1YacDh5!-EK>WdMzkV~NuCw*bI>dMN>y1E&_Kn526HTNgP z<>losUAmMOqe436;}D~2_y3LlFp;jqs~}5+rD%TazctvL^Hmda*Cgic2^4{^X+ zV`gTSlr*daY;D&OMe#i;Bqz^->n(TqkuM51Ob*2gLyXTJELYnv-)~$wm<@eiPVmmf zgqFchaHsO47(eC-l6I)B;_;O_L|I> zNKWHT0xKwwVW~kL?6^}I8meSvm8;BZ^xpe7%oDk-H@4`PV2$rxsXOki?r5Ysd<7@q z>}b3t>xzGa%St(KdlN1&N%4$wLjV<6gtb%ErD%Rc0U={|ox%n@AHx-|uj-H-~RXeVl+NY3^ke{F81%X!h_~iWF-d=k71xd?CvZqDaI!v&D zj8%7y^5p*>+`WFi>+L<@6HM|=zqSOPtxwcB77NWj1Jt-suv4ZpX-tk($E`4niijLR zV5Jl^_XGSnG^FR|_6N$4z`#HfDuNf4CjzINTC5i*HSm#Y?Wsg=dHFNg;elJT!A(u4 zfUWYiv-{PJS~-j;dM&=~e9ub7yrB;_s;a80(QETBWQ&gBj%cdQua$SZpK&{p<4;zK z1ZPc~+ZACz8&dqwFI`d_R>3lOcJ1@?n?6}gWoZbqzB>d&&9iGecc`BK5LVGm7a)B1qk!EZ=J7xu z`TkQZWc}9un?cCm1^M|rj0U%Eomc#C5_j6_5p1!0BROFI<+b%)LVpiQOnarAmei_; zuOQz+T@`RDUpSTJ0ej=hl`G}-r$wd$DP7h;1f{=;!hwAS%Qb77?%r1kuU|0_ZPZ*C zGdcDD-+Swk{j2Es_)k(Vb0{f9-m;`M2oYn}rKe!5J|rxsQ7rnZD`^``d``B(d9c0K zE80j`47e{jJ}pdTF>FMx63o8$)M?S=b}!7UNTR^(K%9y`Qit0dm*m;?3h=E)v5~1D z8kWLuo#YKEp#$R>=3;>4s}q4o@v|hR3x)cjmF8GM{jD+R5ql*Oh86am53q~e32G-5 z?(-~Id$4-oUy8;K(rSYTY}Vx($}ej{hGhti_0owaq%q1+&ZLVHBShSZZNhRG9iRUC z+uy$T9J5H`+V*4$+uXQ~eB&EJn9>o1Aq=|FMix18%`2B`gO3oAMrCYm6WGm48TjyT9Io@>vV0ICqU5W91y@y-0a4X4VA%^6QC3zK6andc zWa;d6QQvkG%m>ocTAGOv`VZ|JZ!mejr}rUvHr6Bh*_+Cz$UnYdIad${ozB%W-r7X$M)7IEbg0K5KFYM>He>g(%^?mKL!MG$&oj29p( zYNcA^#tkS&4WKLtWnxe2ii(bQ{6Y&Urp0wmaY4V{R3NE~E6QOyFs$x>F1pp*@vm${ zs_K7?V!e%u@`>Nx*1UI*W6w$y=U;IjYh`D5aFiQUt%Ihdpx|I z=8eYh1AzUWo}M1Yswyg)nbKlnV#(-sAq~KaFb5D9N|@LwqwhO8Van?gAM$H2jb}qT)mBrOeAB6Jr^Z!>OPoX^5KJf24GQSTQ{LY;_Edi$}mdMCR z0F5QXt$;87Q=G7beO}lG=#ie98bdGjmP0*!u|DC`E4fR7wVxIwkeok;l8zABwqY_c z)C>^ebw_`^EG{!GQG%qpbCe9Toolnp@hW6$Gl^=q>n}Hq%_+nnDMcK zh;Ym6v&@buC@2mO4rZxyzc(#fLlG=ZpWw2>v9cPLx8sqVlESNv8`}9_ZLLC=O!LnT z9w0%P@vwp}yXINjg}~nNOp(5Vl#n2<1U6vqL2AN?T(6$HQ>3v0 zx^_c{4O9SQtIgAO)8!36(yBoHJqPPyx?e!FIRn?whi+`dG^~<@cFT=$-9a6{VV0f- zndN-Jhf}W4%5q~<&ldl#92U|Kf{eMX72j@L!GFTXET7Y|QqxAV2#}RsKV6Ht{)t>V zg=`Gt};U*K@rTTIlOZ!i>7%2gyb|&z6yXZ1ca@v_HFOtsv64*Ue73OB!1SK=0MGd zrfP%H-&geZIb@cq7d6M^R<8T!M0rC}4IlBw&p5k^Vu2d_6UBUmpm+=tgS8eWvfb^R)90C`zDSC5+FDraeeA zzQ|QQV%tcLAkUdmkJcNNzYnIy{AqPSkdttRJ0~7zQdQE44Pt~45w>Xav1L^*&wscv zVu1C|6^$H_Dty+b(}yMZ(yszcH#1gs)J$?^mXtI30;Izi?}($2_m9tEG4;LE@L##K zU8zirjZ}!Tz4zVUSnSIHYJyWUeWct%KKZ`Z3&{ruTUJpo6%`ey1JC^XC_mqu*c=Ad zI#25>J~qM}mP}&>8SC1$YXE>|0*=?Wwzi<6HJHF?E>$*tDfKiR5lK?}&GoGZb?w+1cF*JFmBCfNoP#+rnw?c|~IkFcxY!uInXVDp@ zy_a~s*nyoy4C)KY=nkG^dh@)uZ}pg4of`jVjijXX0srw^N}6RDWQ>AmkOCmW*wB!g zmR9|>UI&;iST95Zk>emSJU(7tSSaZ_8&p(e4;VBf!$d=)H%B4V++V5E(8h!Z-`L_SizK z5gAm!o44l9L=Vqe`DZ$YiCL;-qG%t@F=N`c$JD5d;905|hy-~PimGx`flFa)(@nl0 zUq&r}LSTS9Kai^w9UaXi>wOi~)zfoBOKW53b@bxue9XrSi;tBZ4D9aeqmdv3Lsa{d zCny-Mx@jL?QPCYpeA)hcX?G5xOc1na)<9K)qodZz=-Rc(a=?zT5;bxTd<%DVDso1LG(4-~}F(VFIc zZXTZP-@pHfg;;Ut`EXZ^<@P+Kw}_g}qChAmTv3-fxJ;ZWq% zWJ}9nF2@Qk(x`r`I2&OvH8ybUt%S2Yts*S$BA_agvts29jiAFpwkZ20e{e=5kO_3j zJ6%F~e9Tnl;cz;8!@CH-<(X=wPy~?;Gs~lQFG(;uuWiPZTkv?D9zvJptNC931_6;a3o>Lj){1LX*eg1% zi)ub_Z`|zz0|;u2_wL;zSE)1$Pp1KI3LRbDM;SJh$8{}^%PYh+l?7Y_jNB-s7+4{X zm65LA-k5d8kvH$(pManX#gpza=HX5^ZkmTtM}3EZFU9*5{@})Kj zWsR^BAvPi?>$CH#X@C09@86ZRwPRaqY^a+4{^OlM=KGzL#GHRTse8%1wfeUBU(*CSfSdss& zRPx_FQvZO$tU5ecDI*aq^J&UT+4nd(eFVQ$K=L>(GiwB{RZ2>#56HP69nrKruR$^W zerH|INAtqzK01*K^SPC6eVY*)sl3E5#-nCCR*fl72d1Y3fVfQJ)D#gGHmb7G1dJ`F zckLoGs;0f$`)HhCr^AYs>f$TVF}_Wd;&z;OGLM-A*4Nkf7$W_-i<}jPs{UW4RIlld z05FYbkl6z|@96u-&Q3qGXA{;-n}g0x=ez5PZ=Lc+_->kMl7 zZ)GZYlr~{{Mb&|4e*w7xqX1}621vKew6~@LOkxfi%a@6QWsAP>?y$N$V5)qc$RGYv6wXA}TB)HI*9J zGl(RP`HGYwTdl!3utI@MfFcA;awc!ByD)0ncM^zT9{&^}Ss?pDLa;AgDkyJ0>c(Wu z*R_Rr&e3?NZrHIOI3{GVY_Jg`rpme_kldkb%NrY7sBe(bU_KzN@9eDP#{(gP;nH7O zv;^=PlvTfd`{N-R2FH5x-2?dEu!P8qm-bdhqpecEU%aBPufL!VJnDr=AeKcfY@>Y) zJBg)KNIcKKId|cKYcUacUs#mcF91R}J9KXa5&)lA=iAW8$PQrh^x{2pOpa{g_%eN)LJb0n>G|xd8bls707rmlqR-$gD%9jCd`~r@x zCItWqerSSKs?8f7820jwh|)ZjkP6sjJt=~L9`(t|6f7SH2P-QobShuvU`-8#}wgg+{jsP(}R?T$2ScR!~~n z;unV~Tx?jq&kb)GTfeHRuEub6{^hVbR^cb6z!PDXdH4ZzUW`ix4b>$%VYniTif$AK zVH5|yZUN`QO_95I3&j^hO4xQ{5(N=Na01LtFIA#>czA#bhKkN}V-o1rqO^D$t86Q4 z>v53lFxn%>jr7iU#(M?wbCn%7v)jCe^-X_0841{mODmhGJY~AVwvc zlHa)c+HY@g_2}yA@?^w9-3g3w82a+%Dv++HTLLZaQoOw^XFt`j3+=7)#XNbUbH($) z1EY7Yxi-;yE!9iCAFG}{PMhtwljA+|Jqy5*(@Hp--gtg^VXeUNJh{@?bY9BT>*ro2zGp^od>~0yoqOs} zlg*W$H)k@O!8neB0GPf5!wP5#wy?PK=U zFK`reCC?owa3Gf13K?;w&690+uwccyy1IgKmRjoJ{HjU7?x?{McrJYR-u(C>-PdZZ zk~uGGM59f>nEcomB?vhLL4m2y{QcMI>FJ4yW1zmwd@RCR`kL~Mp(%SAq|!IeTNWm z6riv8Acf+#<-Bt~j})kte}yAyH#lSHUN6~1`R5`?>_f;>(kbUaWFo-Dt++)NQqh%l zqwBe-j|c0iCvV$pD@is$3HHE&zy>q)J=f@^sel%5%X_$3;Jf2O_KT&gy#+r~+Jp@o zu#(d(!$^?Lqc+b^q7xw;m*lFEJ^AWw9eyxE*Ww=Mb&G|B&S)a0i(+Q!!Q5<7bZgtO~V7j$)5 zH4;v_fgoL6%!7gLu;}2JiLZl$11wv3FZ#CM%rW36Zt&oRLmF97ep)&?M8Kwo(j^{x z@ddmk+);NL2^riV8$gs(-tv%MyMA35ESV1;ST&*<+#sz3OPG<7QDycdP6|xZhK8v~ z-8(z*phUfT^$OT_Fn)W2111i*IDfh+Qbonj&})m!<_0%!-ZU_n`1Z~IEjuX}A0JR! zB4T3qpwXn{wR+Y`@>^0Fd3NoqD!Cv_E>oawo6gF1jt1XAL;20K}+9j?;oc} zchBB>A?H~A<_#Yh$vQf~NFmP2xwyU11NI9is@Jb=KyEzRSq4KAIXU?x(7W)pQFDy| zckuA=NJ$||gZqg=#^;gw{d86PD^tBhKJ>O3+oZezIM3m8h9PxmIs z1Tfn!fWp+*Sz1=VsK-5a;LNiL#=OR2~8~1|VvL6$oSwPEIcp6Wi4L8lF^rLV==j z#dFQDqlpLbl-J3@18B|s{`YF#>5&hEgStG)JQhBnxq|A(C?C+PXa8M5M&=mQD>6n| z9f3ng$e>)&GccTEfeV~ku=|^^YhZdORAXshzn=*=6^M_FtOti)K(~6)o|pc0;dk7O zyaaV$zZ#}#fY_mmdf@8%>s!?lfK^e`(v2TZXejDGe*0$h`|G!F^#Gb!6Dv%f?m>w^ z(cs}|ZT;#;QDWjC=+S_}9V{$nz>A$a^}Ra5kPk+14_?A|=U46tOD<5tTx8DYa4`@Y zey%XifBsbcH6(%CWlmveK;J+#1JnKJ(Ic>zg?4qp(w)2F5RY%o!{X0ZSjZm7|j^edYwwJ9I-_;hTpTchBO~fEexCM zS>J`4-Ut#p@>6|NQ+hPI0{oC&Ntg7=Z~0y{2d{-dN|iz0p0wv$YH+XupeJ%|8TVfa zSy^S?CmmE4y{W=KK}lG@2D|ENW)kaWPWa@HcII)n&zhf+2RtPvF6Ld3o8` z*gS6x*}y)8e7HPTIsrQ8ouH$X6esJQ4jd`?tzrpb4FbTkf1zHh8nt(;~8%c9iY}>bAY5>w zM4j@3O&~k)?;~^bX($2cE?i)xj3%R(w6wHz{ry9w(y-UF()bn3&gs$aWpQx;8Ufq$ z*PWc+goop?Sc4i1jzPezlOVeS@@hU@q^so}l?G?597H$`jhKgBQ53uo!J(3PTU7;U zvmV1z>E~A%<~}#9e|S$-wG+xh@Vh|OL?vMQYjN>8R0L?W(7kSQ0GS|kX3Kbiz#fF& z9yPF7pJ@q%?E5K>0r~Qu3O>m45nbkZnO< zMy3(GmOa*Qp`+{^9KyuIP)oqyJ48)!7$UFUO^qiR_M~Tf*tu^1O}(NWa2?%%nVoJr z2EWZWemxm*Vor8+yoUGc=2%|F?1A2{uF9Gk21(bP@^ZoG1^al6RZ(8v;?Ny>5Y@8+ zQYdiG&h@5O0!|LV7RNbo=s}McXMD9hwu>&HCRZal{I9Pf$RMHodsB(ni@X7y= zBd4y)>gwIl)Nfv1L^GTjclHS!PwGz;AgieZ$K=H+$lOBb7{bCI}wnJp|8ldj2woUfv>T@`@)hWpy6bLY+h zzhcu!kC=s+<^7-C3t}T*`V1KxgomIw1^A>gV0RP?7*d)c1euz?Fcq+EWJRW!EvX~c zX_2dPAkWu=zBf5(8ub#;1RRk1R))Gt0LKV{Bdl-)b>RPc#P)w1o|-cPL1K$=q2uu& QI3$Z)SJ6@~Q?v~EAJ#rIKmY&$ literal 0 HcmV?d00001 From 5359dd172724d735041db155326865be7b135a63 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 1 Nov 2023 14:16:46 -0400 Subject: [PATCH 12/18] fix switch test --- src/gms/cp_gm_pb.jl | 4 ++- test/gms/cp_gm.jl | 71 ++++++++++++++++++--------------------------- 2 files changed, 31 insertions(+), 44 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index defe8fc..d1efc9a 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -212,7 +212,9 @@ iterate over event concepts and evaluate predicates to active/deactive push!(active_events, start_event_idx) end - updated_latents = @trace(Gen.Switch(map(clause, events)...)(start_event_idx, pair_idx[start_event_idx], state.bullet_state.latents), :event) + switch = Gen.Switch(map(clause, events)...) + + updated_latents = @trace(switch(start_event_idx, pair_idx[start_event_idx], state.bullet_state.latents), :event) bullet_state = setproperties(state.bullet_state; latents = updated_latents) # death of one or no active event diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 8d84be8..2bce384 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -16,6 +16,14 @@ pprior = PhysPrior((3.0, 10.0), # mass obs_noise = 0.05 t = 120 +fixed_prior_cm = Gen.choicemap() +fixed_prior_cm[:prior => :objects => 1 => :mass] = 2 +fixed_prior_cm[:prior => :objects => 2 => :mass] = 1 +fixed_prior_cm[:prior => :objects => 1 => :friction] = 0.5 +fixed_prior_cm[:prior => :objects => 2 => :friction] = 1.2 +fixed_prior_cm[:prior => :objects => 1 => :restitution] = 0.2 +fixed_prior_cm[:prior => :objects => 2 => :restitution] = 0.2 + function forward_test() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) event_concepts = Type{<:EventRelation}[Collision] @@ -38,15 +46,7 @@ get_x2(trace, t) = get_retval(trace)[t].bullet_state.kinematics[2].position[1] function visualize_active_events() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) event_concepts = Type{<:EventRelation}[Collision] - cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) - - cm = Gen.choicemap() - cm[:prior => :objects => 1 => :mass] = 2 - cm[:prior => :objects => 2 => :mass] = 1 - cm[:prior => :objects => 1 => :friction] = 0.5 - cm[:prior => :objects => 2 => :friction] = 1.2 - cm[:prior => :objects => 1 => :restitution] = 0.2 - cm[:prior => :objects => 2 => :restitution] = 0.2 + cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) num_traces = 50 plt = plot(legend=false, xlim=(0, t), ylim=(1, num_traces+1), yrotation=90, ylabel="Trace", yticks=false, xlabel="Time step") @@ -55,7 +55,7 @@ function visualize_active_events() if i % 10 == 0 @show i end - trace, _ = Gen.generate(cp_model, (t, cp_params), cm); + trace, _ = Gen.generate(cp_model, (t, cp_params), fixed_prior_cm); start = nothing first_x = i==1 ? get_x2(trace, 1) : nothing # only look for collision in first trace @@ -90,9 +90,6 @@ function constrained_test() display(get_choices(trace)) end -# gen regenerate - - function update_test() t = 120 @@ -114,26 +111,23 @@ function update_test() return trace, trace2 end - +# change event start function update_test_2() - t = 120 client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) event_concepts = Type{<:EventRelation}[Collision] cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) - trace, _ = Gen.generate(cp_model, (t, cp_params)) - addr = :prior => :objects => 1 => :mass - cm = Gen.choicemap(addr => trace[addr] + 3) - trace2, _ = Gen.update(trace, cm) - - # compare final positions - t=120 - pos1 = Vector(get_retval(trace)[t].bullet_state.kinematics[1].position) - pos2 = Vector(get_retval(trace2)[t].bullet_state.kinematics[1].position) - @assert pos1 != pos2 + # generate initial trace + trace, _ = Gen.generate(cp_model, (t, cp_params)) - # TODO: more tests that include events + # find first born event + + + cm = copy(fixed_prior_cm) + cm[:kernel => 50 => :events => :start_event_idx] = + trace2, _ = Gen.update(trace, fixed_prior_cm) + trace3, _ = Gen.regenerate(trace2, select(:kernel => 50 => :events => :start_event_idx)) return trace, trace2 end @@ -154,11 +148,11 @@ end v ~ uniform(-1., 1.) end +switch = Gen.Switch(function1, function2) + @gen function switch_model_static() function_idx = @trace(categorical([0.5, 0.5]), :function) - - x = @trace(Gen.Switch(function1, function2)(function_idx), :x) - + x = @trace(switch(function_idx), :x) y = @trace(normal(x, 1.), :y) end @@ -172,24 +166,15 @@ function switch_test_static() trace2, _ = Gen.generate(switch_model_static, (), cm) display(get_choices(trace2)) - # update trace + # update and regenerate trace trace3, _ = Gen.update(trace, cm) - display(get_choices(trace3)) -end - - -@gen function switch_model_unfold() - function_idx = @trace(categorical([0.5, 0.5]), :function) - - x = @trace(Gen.Switch(function1, function2)(function_idx), :x) - - y = @trace(normal(x, 1.), :y) + trace4, _ = Gen.regenerate(trace3, select(:y)) + display(get_choices(trace4)) end - - #forward_test() -visualize_active_events() +#visualize_active_events() #constrained_test() #update_test() +update_test_2() #switch_test_static() \ No newline at end of file From 53f27903f25a1dc30919dea8b15043461b16ba4a Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Tue, 14 Nov 2023 17:34:50 -0500 Subject: [PATCH 13/18] fix regenerating event traces --- Manifest.toml | 4 +-- src/gms/cp_gm_pb.jl | 22 +++++++++------ test/gms/cp_gm.jl | 58 +++++++++++++++++++++++++++++--------- test/gms/plots/events.png | Bin 13052 -> 13144 bytes 4 files changed, 60 insertions(+), 24 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index fbd4586..00fd886 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -738,8 +738,8 @@ uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" version = "2.7.2" [[deps.PhyBullet]] -deps = ["Accessors", "Conda", "Distributions", "DocStringExtensions", "Gen", "Parameters", "PhySMC", "Plots", "PyCall", "Revise", "StaticArrays", "UnicodePlots"] -git-tree-sha1 = "f553bbf8cfdc3a291380ee17f600f818f6cad054" +deps = ["Accessors", "Conda", "Distributions", "DocStringExtensions", "Gen", "Parameters", "PhySMC", "PyCall", "Revise", "StaticArrays", "UnicodePlots"] +git-tree-sha1 = "9fddd996bd0e0fded73a3a15d7f579e1f76cc46f" repo-rev = "master" repo-url = "https://github.com/CNCLgithub/PhyBullet" uuid = "63daae69-5b14-439d-ac6f-096429ca839b" diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index d1efc9a..909d79e 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -118,12 +118,13 @@ function predicate(t::Type{Collision}, a::RigidBodyState, b::RigidBodyState) if norm(Vector(a.linear_vel)-Vector(b.linear_vel)) < 0.01 return 0 end - d = norm(Vector(a.position)-Vector(b.position))-0.175 # l2 distance + + a_dim = a.aabb[2] - a.aabb[1] + b_dim = b.aabb[2] - b.aabb[1] + d = norm(Vector(a.position)-Vector(b.position))-norm((a_dim+b_dim)/2) # l2 distance clamp(exp(-15d), 0., 1.) end -# TODO: surface distances, use pybullet (contact maybe not helpful) - # gen functional collections """ @@ -182,17 +183,22 @@ end # revise for which event, modify the categorical # e.g. another event type for the same pair or another event for one event type +@gen function event_switch(clause, events, start_event_idx, pair, latents) + switch = Gen.Switch(map(clause, events)...) + return switch(start_event_idx, pair, latents) +end """ iterate over event concepts and evaluate predicates to active/deactive """ @gen function event_kernel(state::CPState, params::CPParams) - object_pairs = collect(combinations(state.bullet_state.kinematics, 2)) - pair_idx = repeat(collect(combinations(1:length(state.bullet_state.kinematics), 2)), length(params.event_concepts)) + obj_states = state.bullet_state.kinematics + object_pairs = collect(combinations(obj_states, 2)) + pair_idx = repeat(collect(combinations(1:length(obj_states), 2)), length(params.event_concepts)) pair_idx = [[0,0], pair_idx...] # for no event # map possible events to weight vector for birth decision using the predicates - predicates = [predicate(event_type, a, b) for event_type in params.event_concepts for (a, b) in object_pairs] + predicates = [predicate(event_type, a, b) for event_type in params.event_concepts for (a, b) in object_pairs] # birth of one or no random event weights = copy(predicates) @@ -212,9 +218,9 @@ iterate over event concepts and evaluate predicates to active/deactive push!(active_events, start_event_idx) end - switch = Gen.Switch(map(clause, events)...) + - updated_latents = @trace(switch(start_event_idx, pair_idx[start_event_idx], state.bullet_state.latents), :event) + updated_latents = @trace(event_switch(clause, events, start_event_idx, pair_idx[start_event_idx], state.bullet_state.latents), :event) bullet_state = setproperties(state.bullet_state; latents = updated_latents) # death of one or no active event diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 2bce384..81f47be 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -28,11 +28,10 @@ function forward_test() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) event_concepts = Type{<:EventRelation}[Collision] cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) - addr = :prior => :objects => 1 => :mass - cm = Gen.choicemap(addr => 30) - trace, _ = Gen.generate(cp_model, (t, cp_params), cm); + + trace, _ = Gen.generate(cp_model, (t, cp_params)); println("") - #display(get_choices(trace)) + display(get_choices(trace)) end function add_rectangle!(plt, xstart, xend, y; height=0.8, color=:blue) @@ -78,18 +77,19 @@ function visualize_active_events() savefig(plt, "test/gms/plots/events.png") end -# constrained generation +# constrained generation, event 2 must start at timestep 10 function constrained_test() client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) event_concepts = Type{<:EventRelation}[Collision] cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) addr = 10 => :events => :start_event_idx - cm = Gen.choicemap(addr => 1) + cm = Gen.choicemap(addr => 2) trace, _ = Gen.generate(cp_model, (t, cp_params), cm) display(get_choices(trace)) end +# update priors function update_test() t = 120 @@ -121,13 +121,42 @@ function update_test_2() # generate initial trace trace, _ = Gen.generate(cp_model, (t, cp_params)) - # find first born event - + # find first collision in the trace + start_event_indices = [trace[:kernel=>i=>:events=>:start_event_idx] for i in 1:t] + t1 = findfirst(x -> x == 2, start_event_indices) + + # move first collision five steps earlier + cm = choicemap(fixed_prior_cm) + cm[:kernel => t1 => :events => :start_event_idx] = 1 + cm[:kernel => t1 - 5 => :events => :start_event_idx] = 2 + trace2, ls2, _... = Gen.update(trace, cm) + trace3, delta_s, _... = Gen.regenerate(trace2, select(:kernel => t1 - 5 => :events => :event)) + + @assert delta_s != -Inf + @assert delta_s != NaN + + return trace, trace2 +end + +# redraw latents at same event start +function update_test_3() + + client, a, b = ramp(mass_ratio, obj_frictions, obj_positions) + event_concepts = Type{<:EventRelation}[Collision] + cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) + + # generate initial trace + trace, _ = Gen.generate(cp_model, (t, cp_params)) + + # find first collision in the trace + start_event_indices = [trace[:kernel=>i=>:events=>:start_event_idx] for i in 1:t] + t1 = findfirst(x -> x == 2, start_event_indices) + + # in future maybe gaussian rw + trace2, delta_s, _... = Gen.regenerate(trace, select(:kernel => t1 => :events => :event)) - cm = copy(fixed_prior_cm) - cm[:kernel => 50 => :events => :start_event_idx] = - trace2, _ = Gen.update(trace, fixed_prior_cm) - trace3, _ = Gen.regenerate(trace2, select(:kernel => 50 => :events => :start_event_idx)) + @assert delta_s != -Inf + @assert delta_s != NaN return trace, trace2 end @@ -168,7 +197,7 @@ function switch_test_static() # update and regenerate trace trace3, _ = Gen.update(trace, cm) - trace4, _ = Gen.regenerate(trace3, select(:y)) + trace4, _ = Gen.regenerate(trace3, select(:x)) display(get_choices(trace4)) end @@ -176,5 +205,6 @@ end #visualize_active_events() #constrained_test() #update_test() -update_test_2() +#update_test_2() +update_test_3() #switch_test_static() \ No newline at end of file diff --git a/test/gms/plots/events.png b/test/gms/plots/events.png index 0e70bbb710dbc8ed9c5272101797a2bac552c579..f071e247da5431950ba16bad0bbea8030fe9c691 100644 GIT binary patch literal 13144 zcmaKTbzGF)*7hLM-3U?x0@9!~g22!Mf^;|1A>G|bH;4j)(%oGKBGTQUba#Gx^f~7| z?|IJe8~%_P@44@_*Iw&d*Sgjw{JDY@7CIR^1Oma5k%lQlAa^t%5QJFNJK&RR_3#1k z53-S*6by0;|9RU|kN|;DLu6p0DsGwEv+hoWrj5vZ1;exJPUikWxCyw_cVX09ajGr| zrbI}h6hsoZqE2^LH{4V+1z$h9k2KQo+PHL@5%V=Ai3!7JdTLy;Sa}l?RJyzPl_y%d zr!i*nVKUE#MpZ}8-^^3->yFOawX6?kG3XwvdLwHli^5=Zg}5c;)KKUUeqd}^SWquD z9TfVG6E!?c`X#0q6pByw|NA)K3%B)XNbddn_Yo;rE|n~vEqWzQKWJs*A5Xj3vQ<{b zCgf@&avyrMG+!CkRC{%;k={0stxDoN;(xJ z^G1D^m6M}|LaTSgT$}aE=w!^Q8css@jU_L*UD;YMvD=O*ZR_2!osh>$d$M^@ED5lF zI)$)ltz>b_Z3>6oZKtNYK5%=kM%;dZJ+FvTDsZPhDWpQH!&?`&v0!lTzud4-O0%x6 zL2Att6oOpbifB*Gf+>H@n>ubA8lJ5HRo%;#yg0tO8zfjbRO6Fiz_traem)Je{(U;TlLaTpN-Y4hidSr_Os@ zZH_bNL#|MuFOlQF-`?u&W}Yo(^YD{3&#Y5Iil}tn9ZK3f1JjC^v!c}dm0YkwLtK2k zxF4iJe+)z_>h-2b?!BJITSv!%J%1h%9XyUvG+G8V!&wOVa<@qHU2UQ<$6NPTAvAGh zq(Ov_ciO1V&yynSW|r(aR)*Uhs;$HQc`vkO9V%HKJg}I(8%GMRx^Hr?FZrR)p&6uo zeR#rNYw_z11{SxOx##8Scq0xJVULt1j`AKq5<#5SzYO`~sXQs76r~3_PnQkup<^Q> zdK`BqG0TocHJ}r|wRn-&&8Vq{ha8w0-#N?>HUf_qkMq?{sfaPxocB}rFvoqEwAJt` zqZriNG(Zu1k9qoZsT&B zVTuiv6};(up{UAMYaJeN=fUQ$1_Wvms?bsMptGXO$R;pEf-`OrKi9<|Bf5uf9;m&)A;HgzYqIM2KRiA0760Da%f~i_Hd3%?;=ouZlToH;sgNqo_x-ag z2!6(&%9j#lwTlJorq6CJhekeXor=;ykZIPXWzej75)(gE^T+HKpkT+quPGvW|Gq~4 zsVGc#k@LTi2%dtGeMq@V?`?uV`;j0cMqONn#{QDS{TL!Djo`3&=l<)jJwN;L?wwB$ zkO`6cUKvN95|ec0u6;;GX!5gJ8%9yYQLZe83Q;N+D1^M=#x#$9E$LLf)N>G3Bw4e9o-l@23vT6`Z#+*0F z==V-KnYBSgL6;A61WL3$)O4%vEgYoJwX_;~ZPniv(_QD}-|Dnxvn4hrvJ?dzIF<-f zTf&Z9t%JhCq$U22wZ=htLV*Z+cB>ppb=ql9{@jl>(^L7m7NEQ~U zkGg@yhi*a&=F%P>$EPk^Vg87t|H`VCy*stGSPLYyRphjOz|abj=NWtA(b@84Bv}G9 z%yRu^VPua$Xr+u}m{TiS+h6)LD(zq|sNL~5SHpc67HA6G{(njwaJxPpUb&Vb`b37x zK$YT%SUlmI;6Q16czV5;i%U_z)@^7wJ2UqE6y-FFV*)~dX1K3RK(^uB43o-dyS8hQ zBcs2UEacmqeV~Uf9=8Yli+lEYc`Cl@YQ(}9pXWK>_xq!@-&~)GcxMV%*ilh>A5Gab z>GEW(wtJJ)=b#HoN_Tu6;7<(;^ZLKw1Yq14xIN;vq(jA+JR|EvnN94lpl(qM!pD^E zcS+dQ*j|&Z+I5^T_#3~RCeEF@)Rx7on=5oifr%O6D%;1JdhPjla(hA`3U9in642Fc zQ{__jx(VR;Pq{<(;D46<4lcs`n8;3Re$48L#rO2E6Fbbr5A7fH@o1dxXt_S6qXW(z z?+iZMs`@ICfyI*%f1I&$$22TLtPcji$g6pb0v9b@RQjUEf%aq@!{B-v0n*SgkwV^H zGLqHdr>-~Hhz8yNpNInr7d0KO?JzW}0lq zp-n;a{qt2On~EJ4;;!>UVKQDW00Jevxudm0Zt-V1wI+XhHuf_;y}rYLu^IClNQ zm6P*at8nG#6y$|pS7^mmdAupP^$26J=Bw@6*&VKo{`6XLI`kc~_ix`>qi_6pFZY*g z{|UpuDX+=iFBotv^pM1*6qxkuO+KafIBL&rNo1Keu(05B-k7qfn~u%IBoKKE^Zzt( z4+%xHSN}Gmoa2F2t-WEu?ol{_Y_d)HY9F(mV3wg>+jf?)Zj)_TPzR{xi!RJwv-cPv z?W*e3OpnybA%7QG3rH^f8IW8!X&Ch7|&M zu;(D%+x5*YFTe0C)9dMpJT7%Vq%>dvwnp9fxQstlV~$TKYU|r{j;0{*^yYyvf$#jN zzuY0)P3J(E_lf-B;ZtE{_K9?1LDZ{waGgLxf|2p5iVC;uyn_31o;>_jt5hQ2K)SAt;@lKU4< z(8SK#7~{qOZ%8cFlaQ2LU9zcrr(NekvR+3q5tU)rOCU~UT|y})USTM{YBXKDHu)k$ zAq>_FCqE9flwLpPH6xqN;K;C;Jg3WG-9!roZBYscNE4? zN1q>dRj0cRuq^W!oZ2?GN@XBh8a__X&20e{hlm5*6&bNv&xVKf(&y{m*Qv+PF(M`4 z(9+St^&eY)LRl$qV8bs!z;iA9!8z zFKG-!0N=fyRh=^0M#kwqS|7Ka-gc>)EVdAY<$#LydTTavKsd1L-Ox`h1XsjRY{PY$ zq*{$8D&JyYm@<}YbbRZcHZI=BsNWp+U!KZi=aQ1ek89-92yvnlru6ej;7FbIBR6|w zeZj{u`7ZD;$PwUjVw9C-A|YY#EyzZ;f5oeF_d6{Q@kbei8QEFqtq*^2G2W*{H}cyz z$QTs@q}~lLjRY_=kj4A+w*rX2I=5FsWSKDv0g- z9orVw->6JYM=@D}n847rtYq=3waBwchmlRqB62+Zb??nJD#ozcCTyh(T*U853?iHP z1-NKCF=;VgOS}2|Xi@p}y&xDFAnijbxN&pw$nZ7c!JDwVpMDd>;9||HpMf~Sd#s7z zm6BpLjfWIIG!I={QSJ0cjH zLI45~i4d^OzR9Zc>CI~vUP$}0zpn44aaX7PITuP2{zsbb=cae-246lViv~R;84nD$ zkVXV!y6E{Z*y4HWA@pWt#UwhN1`1WyUghl^E>IkJ^Q6gq;LY|_&2TLP6bcDTcgs3(4`L$1!9Arc3(Ay!BxZ z&=h3gSy$6(`-@i$%6+w@zuy&#gG0{qZ8770G@v1nHq_2IS}9@@5-MRYjUwAH&<&hK zQdeu0fHwB++o%8f#5lB$XQk zKBu)GQTHARxbDnU7)xQ1=k$JB?d%BTp1b_Y{Cz@0uEsSfk@{^D1M=40+PAKD-^lUX zOCCpkTE54FnaP>;uPGm5i77sUH(cNQ};F_sIldqudcts5F)%fK;Zhk7?VnUaTq<^=I%(NlIFd;^N`Or9T59k!0H7yg8<+si{!#>RT!& zjh60t9NEWMTsNVYQ`HxqR%s|xCSq=G)F^#?e+oeEA+P14)0oRqVi=(bnV?H?->S#u zb{z&f`cBK~Txx3S-kyC)Q;wwRmWI~>TJzk-ZgdzUigoLBPme}h>p|(Xoz&BTRhenT zt*7IqZfD5=0p1DK&7o5^Q!D#_G2r0KQh;!Ba_+XA>TQmd=zFZc&LA2B=(&^anw>ytNdCmcL;{22 z1u;-Y#Kgs8CPf`vaC&9Q&F=&~~rJ81^OT%@cStHCU9=X|DE{ z)Y0-&GpE{kVNMwOP;mDw}cSFYJt15bdR z`Ang-#4H*HHCo(Tpof;`&c`C>(W%?O@jU|^0br&f6C>k~Mik(He!w)5up%WE1OC!L ze6tZ#FnGvu@bN#2`uak`Zk9L8Dx9rWmyXua^?e@MP?eL6@Go-p4LNVwY z5Zlye<>{a{X-!z4jdFTGmC4~9(-VsI^`ywQONxWTq~iID$M$zvwT+$%0&DvyTvg^M zIqw|)vKm8+EB_;k++L)B1^@w~kS*D35zNf{OcWt7RvHm-x&I;}f0Uy@Q${{g%(Yo^ z^-@O7NBOvnPn{l5^)hPKc3c1*1rbWaiTdzmP|(xK%>^};;V-Ms?(0fWc@ETve4n{o zww3GicQ;K9X%q;VbQ%m>FBN{*MdQDOpZ{WKpq2*&cK8YIwD&~Hjy?NJ=wNsVp4A)k z?`HwNASMy&`t#qyMr{W+MQ2f&UK0hr(e)aw#S%=O9EuNi9jF_$3{s2fw)2X-UuV(P zPxLiN9bX+gc+e^te*!jm!8Hl;S;PW`O??XRfd*w~P$7#3CuH@=9->y7TZ)4^qb@!acJVho+q4G`KYUH@<= zoHO_hYPhc^BF9`gHCWzr)MTL1Ymbb@O?`Y3`aHgRk_}H4=n^$8XHc$3%p_f16i+pD z-?+im6017xy50oEi|H?Zy6rqiJJIfeSFEEPYRu^9g>CWmU53&MC%o(7AK1?LX=|{^ zLkzo^NaSQ^oYBNIC5zNhi(VLwCG|&7`xUb^dq4kDd8~Y#8C$4RT7G-2zndFVIE@t- z3zzWd&`eEbiJcH2Gq#EAb|+{>H3ygkP;^Gn?frc-fb1NDr~HG16c_KEiD)E-J@heQ z95HAzI=6q$G{cw-fvT1z5d1=wE%bV$h>=~r(&u)&-3*vr5CSAGJlndS0FbMfD3 zGoHym(a7!hYbY&~6x+mXsUOjQ3KrO)^HZ&)nRhntkERjt(=dzfi`4vj_;fhj$Q)8P z1Qy?{SCT=EkLcQr`#0SZ9=@MYW!hMXHA*-#5Y^_<_zWwkWqk`2`Lddqc;aNwlxR)2 zrh*Ywd-Zy|4Tk7^ty)hh9UYlSc&(pWd~D+`+NdL#2&jJi?Ou+_B(j{0=;_Y5tzg#o z?^)4*nf!2h(!9cxH%bx`Yw9M1e%H`TXlc1@Hre!O1F2|qK3z3!mJ|W3B>EckbEGn2 z77^*=@rLuEmCO9X$s9u7|> z)!bwY3u>o5$a|pKYaA0yL_Bupr<8lA{`lWA2#P&37Ol+CC#q1dgD*8bsia+ zn8H_-P@#ATG(5b_&cVH zwr`F$BIfi)2AWE@NxtqLEdr9t5;Qh2=uI-i5rYfvvE^m7S|| zu!W8l;wXfKjGDJ_g@tI*JO*fwM$(F>3aRZ@QTRyr|7ght*@UH<$Ph+$DW8VUq-*mE zQb$L(Cyr3=-(vf|7aDYQZ&>p0HhHWJE}v24^oIYbGj2<0It@7{{kEb0R6tL;YGg6fcmj7$TB70JHja;>D|w6#u()?o$( zP#Fz{D33NkMwGEmZQH^-T`72{tNddRRFfPN(fL-@6pROKaA{04?RTFF0{tHnB;jLB z+Og+ji6rr)@_hDZ^aZq^ z@Dp+$6+?WO0}jiJ zeuE!C&uj$m>n%E}g$L~+K_??WPX^t~ZeTJ^c zsf5weFRT<<&hfbtgb24X*qFb^7tqQmY3V4?u!4uf!%%B@3JN9}EY0h>tj8MBawH=} zr3D;+Wt~<0M97p&$tiyJrT=&q$y$8iP3jZrFO|;b6$>3GMcLk87VT5}&6a(?Kd9Cq z6b_9KV*-MSl$cm-yz;efWfT&I0S99P|08?yaw`(O=@OI4#z$nFu53BE@T#8DZ>CxE zNCt>SK!_tSRiydAm_M3^*6e~AM+rp4{~|%~;%LoNBKSO=rOOB|j4nC=p>e#@RHLcZ zasRs@QR)0P?j9o0RCXI5>NM<7dRHHuycKQJpo2OCquPWoP=E7JD+4b3HQZzh^aZtz zoi!z7nmFWotHl2~*o9A*-sclKH5QQsw18(`mAMos zj&v`gdY`W{`<{sxk43LV4zJHtTMSdU&b)J`ZLBi>^xEh8LPA2qVm7Q*r@?V`Z7>%R z1-tm#@GHt<60Yj>*?B0WljP8s7ea8EhiuTB(WkBfgSpvO0VcGqW?{A)L+Vj1>&!6K= zdjajMq?4hQjQX%kfPwBIkTsyig6v!ShFxMip;4cUr-l^6_i}KJ=`9Kw1cNq{cIA5;LlR(5cx znBKM>D^3Gue_tOin69?SbtY*{$j-K_(So0*u!ZTHXx^Ov%1%f~sIi5)e$N*1?ZV;( zpa^bCQZQ^$H5r0Up>^4nhyy0HIn`^Lt+&5t2?8}{X?fW@vXaGaw$b_UXFpE{ zE05X_3~Ne$M-!k-6)0xg&NVmmn;n*1trujgmFj4!sEjl^TLpHM4NC-!oNQ0mSx-8i zAFQsdtVF24ZEp5j{hBm6If;gbMpUMntFES&oRpN5ot+&U``r0@uEle_OrKpw*~rMK z=j+hWkX`Hf&vt)ApfnizD2(LTRkBct_y`FKUWSm{qGDkkmpC=(OxN06?lkS9qobEz zZf(uA2-4CbK3(l&Wd5}rBLb{D|CaVh&}ogsxEB+N86~UX)vx@#ycySd@49&(&nmNi zVB;_CHk;{MpgCC#=O1rRyRQ9sM;RkEEnHLSbL~+zDP-6g6v6$XKw*32<8W?_@Dz}7 z6h*F8y5_16Ix*P0KYf}9w{x9!9tWz(FYR;SiQe8`;fuB0j~_qc2%o*Jr6?TOJf3qN z*VS!wN)vRo%@+1X#w1;E5-`WLH=7dJLC zl8ho%a=y;Kz512Rsx5pnVO*kB<8eHy?#+C7cnB6R=(07@wA(7k$+`b4+jlI-l$(%{ z@Ug{E_v!p$jL3~fh46z1@Acb!aGo4JFHp$3INjCP)7#O$1(wq>Fbrh!S4l(?M}qls zi3$tq6G}y(HK~{BvGzIo`1qKZtPJK#B{5Mt0om$-O#JSNam*7Qiy^xCZ5r_q7R^dn zvT10S!*VolsYZoZ6>-fAEv=~%?K;gWGu69`{)i~)k1ZTmzY=LCs~2m;#m6h}Z=e}1 z_eKNxu+44XXf!FF4#;v|=c_4|EF{2co2etP-DKP!q6rujf`Wq1=18$aY3b?Bz&a~g z+V*?N#X~UnJz6SRM9xy=b>n%ye$@>WW~z?Ccr2?>cOV|#D03xo>4Bb|^? zb9wpmSNOE6DUg*QSCKserC{!~y<;JPugemEe4 z5}S}fKtiJKP=bpv&*9+UAeX{M!D-YD(EjpdyUBjBb7VvX0`9}ARlWBl9|Ud+yS@wG z{dAYYw8Cr(754<4Q+$9mpL6sCO@T{(4P?Us!_udLuxacAos4RXBw4g%%%@zHeQH%L zEdh|&VB!)(&WnBP9?qoq06{+qa(`B+SG#@GRO^l|L(IACk4hM1oNw`Rvrz67@Z*Sudci_ zlf5&(^!06y7Nr39D!IF#*E=i=xNOy%e8HbjKRa_P8}|gjT3^ox0l_+sb!~anda+Ri zP!m)_jHfOxE?`H76-8QJ4(5P6M)bw?ReG%ty)O)$hsCC*l2TFGaiu{pDFluVR{L7( zKCSkpwCI)XvTNz-&B?=cnSw&C21iry4uG3E2PwULdbR8B+U)l@J2l=LKIIc;sD#(EC(EZuUoOgSTg&L?mC|U)O#y-1mN>bP#Cd z2op`N?2yo|3V^)x&EG9(e>w^5wqEXl;m7-e0`=xw5bRCPo3(~r$dF1F5QkwW__L0E ztYs}HQy{%Vaj51>b(?tW z%S(1SRYYV=KF1YlTU$22=H_OQvFB#SZ95OhxJ<%ATA1OPtZd-LjKYtm}GG|6=Gltlj-n2MQrcw7LJ zcQ0q?>FL=t$|3Cl#Ba`)V+OKCZo!b?^S5s^ZN6=zmZOcDHqBsv)6&xPy^oXZ9Du26 zMb0aqqz^2=O5)A+B>_G@n_lw-z;e&s7G5467m(;+ zjaMG}mo%8NJ{|O*}vR1&SU(07U=?p!~`uvs6r4vuf9Z=^N24T{PvB zi`RlXkiB+9oJ>B4r8k~C=f6a5`Ry0l!Hoza288$f!7t!>0M|f;wwWogAKRI>YkT5; zVA5ms?Afy()KWk~%Czg)=;`6$PwalcF<&|up-b%2bPb9*d`!n$H^XK2jnu|-7$z#( zAw}h_qN4H`JOZpwy;#4+;}~#Cty=5zAL2M@XLc`SWWIou0V_{PO?3k-R7J;n@#_3A z+v^~Fc)fYEw25p#`>YT?=m8iW(L-+6Jb0?E7oP_*j7q@i1HiP?IZwNucO77&3f%N8 z9*l1D!O&1NBVTp@3q$r<2=V2!EP+v@issv!OVyaDleNJv*Xr-f`z+eETL7H^%ehQ^ z_`=I`iw}zcbh7d@a$aDF1bQNir-N92+*;0)B- zVEWGWQY5?C_sq#k(?0grGqL{C(kc)OF`t&3ng#JjBaKLA=`D-tpQS!`csB=Vq5w1V zrti%!UnVvhauO0=tFht_9~igE?)tS|Z`%PLj476byud4#;{hV0ad$|5yq`|BBqt{) zHZ~Rvi6gkxOLZ)U^1|(s!#5KF0ZdHn?(J3cv{P+G6X@*fDpD9005P1-CTUqW9T^dk z4$6mfX;ERRUJLqL5F~Bq0ryg=Sn4JUu77dV`Qm7tWqNJn*Xgw1< z`+TLT9HTlhIePA^Nrj3cU~7N869@0{@$rp}jY(KEo)2$%V22v_M%@DpS6EmW zBnS=7qW|5xf9xY}(ytP{<^GWENHC=;U-X)9{Mf?6FRTk;JoS z-!6_f!Q;p&;IAiMXZC~iHxCXLOO>F$x~*OfOO6w4DMjVwdvQ{fX^+gm=SW0?5f;DQ zT=Tt0W+rB400G{D69ZH}XN!^St$Z(k0EU5#00O@o`~n8Tl>j+GWj@{<7rDJUL=IeB zP0_zO4!Joh5CPqWIryHTmR7nQbpFSW=a4c`4?wct$NmHW3IXtfyH9OxcL2cxlpiqt zO>hF->vYx?v_MhA>sOcaw|axWz{v(akP)XKhXKk0{&+S~Vf^UPBi}>jhf2n@1WH*? z4nEUDKwWIC1T+zl-jv-FP!7=yIs%Z0PUirKK!EeSA{7cA0p2#(?B3dQx!Wef$yrff zU%#|JAbe>JIK`*s^W!lo(2xM^up}o37z)-`SenLbrw4M+*4DNs9TfMp{-2$lVw%Y! zms^h@0GA5u$Lq{`B8b|~p}5>!SK~3#?yY-=7W1vwpbjvbNn&AP*$n6R1B&kzNzP-v zHq&spGuy;WN%fn-9RPw+ke8>UUJK3&0K)qQAT-D#uwlP01vC1&8ve z1PucN9M%L^mZY}WEp>+z#bk3B1`NLZ_We7hh>us7k5vve(2|w)M$_n7iRDa8Oh&Jx z-@SX+YV536#w#M?3jzkrj9Yis&(F`>usc)qJY?eY;111loePhO$^`oXdaVqwVjG+k zHeDt*wi{=F5bWwKy`Mg*&+ugQS5#C00#kNsCr7vf#E5jzG=U*&JpdSZZ_y|;XI2G@ zp%n%Rn-GX-K>A47bW&4N>Od}7Qr#yaaz5>OXcdLVNF5Ez69ASqt1@ALC1`RokB14wsaXJZ=oH)v^o!e*Rpy|Mf;M##K@1-k%u(gQR_e6*c#Dbxo!8rjM>~?3m0s}; zyi9Eevg_1y0I&gsB;m^lIM)IWLxAv}cO6eq6y^ub6T=&~96)S9H-Y%0)@B;?Cgjn* zpuGX-9_Z-kq9_DUfNNDT*fc9w(#+q09atMKh#7nd1P#+3cWT(uLXv6NhG?>B*|^nX zvZB(XtXsVltNxKC(077;wFD&4(Xq0}=;OK}07j6uufQpZBDE6R`o)k&6=C0-t6&Tg zC@n4M@tGZWK#XP!xgP*l4Bq2V30rfJf(o_utL!OZ9uVk0P+}r?IX^^gqy51EXfQwj zIWHA4CFR`57X==w3ZV4M$`-Ecfpvny?LkCmX>x``{YUDg%c-ue2273@XgiydVb$OUE=|*s90qJh3p-VahgrTGbNhPENM7mo91VtE-5D-B+rSrby zIrsa{Ip24Gcl(FG(LH;=?|RpI*0Y|qCiUu9wV1ct*s-^Q}~+v{Y6&%Mt+q*HPMNQle9{UcthcEB^7)* zCMGRGOh#g24w~>pIGikARtXA?dC$Pi#3XM`L_(aDe-}F%t~5#&0fiRY|8e^y*7sPx zb^}Ahwf(U(aZNI3;jD38-`lc8%mF*9>4{mfzgAedUb+is#rM#Qvn&(FQU0h(e=oQq zuGMs5wK3cDOW!;Q(d212nL9CQ$W0+jgpV0M{?s=oCMITK!Qwb)-oCfq+3n0YzkelG zd*X7fcwWlKcbaRP>!mmPE87o5-VfsR*5z8CwS-s%>e$%WJbai4-hV4=m=!xry|`u1 zvh}0V&v!L&?0ds} z_9p+-wQIL4>!*ov{Hv>xjy{s4q3Rm^f7(S@E?KOC!d9TW8+;=JLW++W;b(3)w+^%S z^Cs~_=Mzc}-)Zm3l^XGxa`HW~Txt0d3~9fwf@@3QW$sY%sLm)hW?c2Ts9t#ZT~|6C zPGv`?WDZQmfiCmwm`$e=W8Y#KlTLdEK0+>nrs6lkO+Tl}&>!T+m9Uq_FeVFt8=tUu zTRH5NzkxtxL#H|!QM33wR0(^ca=9P^>e9u)Ol{Ock^l9;vH=yMo~JP3 zAVN7#Gbw1i>c!5tvPlYV5td`3e>tjTijoxg4^O85oXx7i1iSV4IV1J@Bw=5l8a>4x z?ZcS2cf6wALezMeWAarnI!7G^*|DEINv41@o8!wRzuM`jG9DS~sdXSiXN#K3Q+rry zBpH%|^|b@lEXW^KS9cXLbehCTezix4Prv=@4}2}; z$|DXKA&@q@{%2reMxV~FNLW3D{A2+4@hB|% zaQi)8VL)PH-A_bta2-|fwkP7KG#hT-jQi0~d~@?G489a2t2F4b&sl8 zbf@MfPb|eTU~{Ia1bJTw*3nG=G<{=f3EyoW5QMijU-Tw0TDcdx8U{K(3Q$XI%Lu-vJ_s?dcG0{HSVlpNfgGK!lwrCkS zo*8XX1o!1`kS1MvkGEvQ&br>oE`~MHk*HKZC*KZ2qh|*EX9hDdTl%Fh4e><8#7(4sMN(loE%7$Gu3#-n_gk?~wKqxpfS0E(ZyMf| z**4e?Br$ywe7W0eq3r1K-l`fZU+S~z*4|*P=(V#w1qIzwIR_BlH1tBkMEU3Uz zC-;kZH^jqUZjBEJ!&1k+i8vr&^YUbt?PqTzE~(LDF?8pS1;aXIpqc=eQlI45fTZ-d zTvJb+{0<$nd!F8i2Xcad*z)BBbvs3NF?bV02a)(K5klwMj*~T9O+yaeYp2lVozz%#_>hmz-S3Tn8BnM}L7iJI^N;*&q# zpffq%L=P2meCv{9-E^+JK-K)XPwE>B$CiSVuxyc5`MO%4pr6OLIpW*9CqGDZct$C2WV-Y0TWlCZ6r z9KxqYrp{m1yo`9N9x^rjuu)iY+h@kbwGi~*zU3`W74Jbsn*E`;IzF-FSyV@8soy2B zABq|LuNg=4mc~&kN+mUV8oaO@BqM@v&5SI41%^~KG51^pJ;Id?sUn!MK1qfj&pGau z=azJPb~;vB?sb@&YYDJ0i_Y`LM4oADpkwC0grK0%mBlD&BJYnPENE=(t9mEHuBrQU zkyD>?33XCZ&D`u+8RwsY{{gFn=GLP~_%wj6zp|E7h-_jrD$}Wz9d&XzwEz*t#G?Ko zCxn1WjC#z!>DB9cVszOGy;7ZLGP>bmL^Pbv^3ZpE49sE~(znKnCo3W&^IKcj^t~}( zLohH5^cUK{7^#g}T+Scn?omnKoM{NCrD3A}QK6aq3Y#s=s6o!RlawnbYnKhi8tMG`H^D%PfOGik4O4KUhkM z-_8`iVIGK*Ljakcw8*~PizcPxsyMA|5-u+T(5Q$LJsXA$dBN^SVFTz8TgW}N8=Nra zic)HPOapE$N?B%Yr`vyn>HiIRt*%dmAf8rv?%nVbM_(`EE$6TmSoz8F2lA=%S&p;! zdra#s;1w2G9yucc90AiUN5JZSbaf3AZhX?6aEOCjfx%rt8r#;Vs#TwwsDJv*;eQ_(oOzFDF4M-)-EoHe)4$} zPZG;l2iv`f!7Vw}=Sm8=1Jc;wU}D;GpngJ+63W?8-!;C=M$5tW)*^BB@K2fwz4B4} z&2oU=dB4$nU{>G!d;9rKjn5*xIybaJFVCwbJJeh)>lU{f?pn%aseifWU8hz|n)6%@ zTG_V|_eh=W)nq~I>;Uh1HmQ(MiPOv`>60f!eZo0U;LN?YSbt5#$YkRj!R+ZTRaQ<+ zw&bLL;~`obovd?KQ2uL7|Ej!>FgQ&Q047>762%1ZdGaF-3y!mcu8X+(_94Cwi-F|& z+(PMC`SQ_{tkcbf2%Oc4;j|s~LUzC{#Oi2Fn+4zTp+b(fdghM<-oE_#Md+7kJm!-p zz6Y8O6!)3!KsbwQMhvIZ2QQc!<_iM{|L59dm6m26utr2EnW880S=&+nhL*@1&z7sB z?&2sZcg5)-Q;jWYZXij}nSCK|gwhrABpoL#ekk;NC|Z!Xo!pU@iHZ(&*f@9{6KCc7lP z$$g*S>wK%xt8}f0sS3h^O{mzQbUZmOm(O0N1JOA`z_HOWl#Pgx%uqS9{i@kTj1B9@ z`!1NbM+OKeAt81?*RdIgx|7}e_;i~Fx9&iSH5z_O{^ z#XZ`dG}qL85`)Xu5_fb>tNF-j%orVg;p2;D#6@7naD|wLx_Z_iL;i#%E#p9*%JR$# z8&iYTUka>o<9a8*<8MA&T;$~By1KewD4H&j|i42X$7tkfWY0Ux!3{1?COd-5Z zRlc@^VUhXGnFfvZGPnt=81;)6BQ@)ehJ1VTDDAin$P(kzD6tzsMVo z+SJf-@#*0#5sh5Qsdj$PpC6Zp(kaQw`+9qoj46`JO`0RVwOH1OO$l}4LOL$Oq@3xz z&g*BU{V^c{M+;9WZw_&$UAOGuYQJt=sD*l2PoJEeOs5&UG+}g>Y8L~*re$GySZ?Wc z6dHj6v8cmoS;B$owpGs4#48kx>gegI7!V~5oi6^k(QY7+l+gTo za+Of#!lKTIl$4Z(GZGH}Ft*-KoGTsVe|2#hbat>dQISSC@J4Y0)?fagza{VC1#gE%c`h`9{;vZb`rdY0^{nH2$3I87)wZSFrZ?Z`MH&HJRv=qoU<#)IRDsXUD!1 zF#q|^)m2I5PZ99g#@Cbfva(%?jQ1ZL8BW6if-Gb*+>{RV1zV7neFDA>cJ@_Bc6Dc` zYHUqKRkf%LX@j0hhugCqAsu~Imi?N^q-nT9LU*XLH@_I6aMc=7WqTJTJs*DGWdoNi z+@|dO-Vy`gEm^=o|)o-?%^6S&!%eQy5V+LJ5!0H}v zmJAL%Pjo3O?;iVj%0l;xw!T(dTUl9YYL4#x=)ZmYcFntYIv4j-A2IWj5XzV|h&zzf zJ3Rp$6jLA)rK!nxh&HQF==w~*rV?61J(O;}ZQ?U|pV@eNe0&^_l9!p8`9Y4j_gt&r z)^uG8FuVwVmZrpM?b`Zt{u(a~p|Yls1RDICk&z{8wY0)ciOrZyOiTkMS<(@J{1EfI5#8fdQPikpB@D3h|=Xk@EsK=%J55Beq92cpkLmt`LWexw#$e=1M2bav-L1ympjPt<2fR;Yk#}Xr~Lhl@*Z$A7&fQ^yy`;=YYw`Omy zD1+NfOv|59S3ke%nZG&`Ddl2Jj`We1w&a)H2=cTd8Xy5WMbn~H>v8U{XKZN43}_%P ze0)O3SI!Y;!NaM=XSf*X>MtHiaOebF0}!JZ=Vy-GMoW^|*X#^vABK8j?t^Z_kNuGu& zRzCdgozhsqwerp|?u~m$`7eBr%Gptofo$I2nae=T*9ZxP|E4PO!_lu1R$Ztoe1cTF z=^dg!UohqW{Cu=IRnzqvSlYwNEIFP&xN{wG=SB7=F4lWgEP@4R*FdsYX(5H{UcMLZ zz3-drshGsZXz*`fAr^=nnFajYeleu71XD7@IGokhpN$KFs8xRNhlkYDEr z&A+OK@^+aOFnf-UzMWd|f@9~9zr9!fv?FdBadAmZgsfYnp8|d-Vo$`$wV|h9+1pPq z@mwMYQ~l+Atl3eud{+$(;b6d z-dD6)LCGjZKsc*J3lZ-_MHbO{&w-06nFm?_{`&Oebfv128}kYsVkX*Jsj__XX~Z5{ zWfo6&nh;`f;+Bx6z$^mzNHdqd+hWYPyy)nx+r)!EPn7h_5$LiKDJ{j!UEquNE zMh5ZvS7XF*q{I$X@40ndbF)J6m_uEw-TkGtGerX+!O@4P@SPe zDNw&SNUrDp3tMk5dyP3=e|b^n9C`yL)cQ&_hCQG9wi9X6#$D_IHk(HaMd|B$!?slM zJkkTV=FtK^SJ6}0Hdek9UL8d%t~@dqdZ8npqf+3or|c!}wYo6f41e@YR!O?lRN`VU z!MprNzxB{_iV#`+w-)F?xose5suxecQ@d@0K(KSlJt}nkGqj%!px02_dJdWCqCwso zSR9;9v1RfxKq#m$J>im_c+(0XtrY^^`!9t7)Qd~^3Z8U?#7Yby3-{GS>?2PNLoJXK zt*M|q;8rn%cjy=xJi)pNpPqPA3E}3vMhWc;D^N5wGvWfWaJ+)H#Ce2XES*&5p`46N zZ*ODu9i~%DqGV9tGjfzJ^GOgur2Jp@Gef@C;#9Tcw$0TQ)n!tYoZH$7=9gd-6O&AU zj09W&E9*vZFgsUDGeMN&wkr)jK;5zkrJr97L6ZHu-c8Gx6s~}lg2TwDHTj)k!j<*S z8z>j%0&D(H3!-EVPL8hPy!u0C2eDPWrDaAt3`8o&|dIXmlVgi8z@INImU_00+nbs2&< z0GERb2OdqVd1Wfl(Jft&KzVx^nDYv?gt*ipODopHG#L-nAIqct2)zirkeIJS2&q&l z8qr`9x}1GjhiMj*8jJTZ5!aytz4@xu*0zX3+2qTF~&_7Lz>2g=9-cW8hiTH8-fS!JY?ux3%2YWISL(+p;~>~ubS36$xB?;I%!Zn|;L-+(lF;b8G91I4Y< zq5kuRrr`is4_-Yn$tgsrN{6Z8`+rm@lYdkwgC3KnjlMt@xMFgnb9Zm7YXc<2x@Jm; zQ`fm_G?54$OV#)fB{woT+3nP0Eb3AZ17c$%y77~m%}AC^_b`tMA3uUY$$%g!8(GIZ zqkGCO9F#!}d8qr;%JqeNe!F9C5g(r-;HRkn&J#4WRzNoaNbvW%=_r-agB~e}`)}SQ zqZA1;*g;)TWM;GWcSY7sVGJcYoQH^4h_Jt?Ce%i4G|S6s^v|&$VB+HL=AkyEk9%Zd zZgA<1^AOTTzTr;LyU%RhDx|`%qsi!v;WClE(*hlT@ZSVAgH_g6Hri>E*eqHhB^^a@!QVPF|Dar3M{z;>3GWHr| zuD_ZD#2i4fRQaN?E=UCh-CO1@ zWofDJs&GQugo$3cZiU|Dh%_&LI^%XXb!`WLUQUw&CoLv6ct}5*ly!Xwg zbKm1*Uc#;$&IThg4G$y^Pppb`pOv9NSV=20W=$SC{(llm&03Fm722G9mDx%7v6M6;c9@Z=$T*MS7 zqoehwC}@QoaZ16Uu)800F1*ut2}#UL(3NSb$<%`#fWFkM2cCzglp+3auUiYWI*P{p zo!OI#g(dzAR2BV-o?=@kMMnqJ8$Brw5GKRNJMWL@C7xq#%VK^mCj4oKdHx zbcsFwy#IxHpX@)SZluPMlmz5lNdy8LXRrOWDwCL1tfb#G7Ci%3xr-1}_fEWse=%xn90-5@r-L{H;&RD&^Mu3q9T{GuACX%*F4nbF6kRux=u~ zLoFIe*U5f6Evq$-{vYwc)KE z9arnOPxy{bYWJJZv*l3Em$*RA#ohgSE|RXo^1lDq9uAWR$U? zo8Qv8yOLS-^(oLoVILo?)X;d#VBEZ) z&fD11(mw=?G2)r^FGu!V(15U18Q|kIH~XQzuQShWovc1*-(Gq6m6n6;2-(?jIXYoh4?<=&Qe85UdutLhUgDM8xYV@`xkuu1B2WDe6RSp4- zTCBNEosodRlwsrtQe;?+juQ2~COj_)l|ef_HSqgtn)cP{!$+u|UaYx+sq^TXmrDs8 zPGd!sjQ?VI7+3E6n<8l1=JYhCrnl@ZLB}8eR%hWc1+EYXNfB}u;BOsw07Uzb?g$3O zP2d{5LUb?GZ*1D>r1}vwuZiyz%vJQOps0TsCzn3u9O4otz@1T#0C=2T}B#>TRWJ zE$Xl5<2*b8!C(G(IB3x1VFQD$U!49b;_9c(k@VXKrLPT<*$%Rp%ni9ZZP{*C1tS7O z4fDY*{)g+V;fcePbWFl4O~U#9Cj$4!6^2{ke;sK^8*=Ti7tB03IQabeA#yI^5{M0~ zDWw6&+gg;1rch1b(O7z7(%Zl8DDW#Oa{ca+UqXO)Uaeo1zEv)Bd7$m)<|Zc>R;%krTQKUoyU_In4NJuo zJ~=Qvy!E9>)z`NXOhTkDE?PlOPt`66+=GaCjqABtSmcMJiYs-^%v`Dsm~P)L;_cnY zmJZ4cx?QDjVzRljLr+G=?9DhQPu5ps_lcF275PXppKX(t$6KNg1JG=1#26SXz}tS? zA0)V?`zLMj|D?aa)kp_3|1}oM-0crZx3h$uIr;eP8Y=Hc5!l$-7rQ-2TQeB4@G22lfQ`*^7 zCvIs~vc4bv2`kIXb&iwjRg*I!6XWBc{$JUs>fq(!adz7R>M5Uu1ns?V{p`&nSoDEJ zw{A(EyS5WIzw+u(BtAVb@YJx;>0`38D= z&$F{715b9>#)|uUdv8iSez2}0WjXcE#=lLPiOK5gTg$_Sj{QNdTjE|B5MuEl8 zxU^K1l?}FEp04C+Ja~{TV5bgtn9gHtJHPslq1LDZyyR=EUk!UqUtb?kVignUyb%LI zBR~vBmwn6k3BbZ)Z}idf$B#&%IX2&(|AkBQ+UvuONmadkFHSBl5L;5NbFI@1lD@kZ z_V)J5%F*mG8`p?? z9#6`qUV?y~6Gjd^ih;ww*4QPm$5gc)8QI#t9=0Ww->dGY%Mx;Y-uu+Qw=b5Gi<5KC z_ecDhtbvgc4H419hYuNI7J+l3qM=opG%>|2f)d};W6#sm(++?M zaN^q#sJy(*qeo72ty8*e$6|-B0q^_z5aZ*Jk~eQ=ff5ygLkc-LJ&mUpupN7?rlrS` zWS8q_3;ge}|G`wP!x-2m5c43jo}Z+ECvvj0Csy9Wd1KtLeRaOmV89J*#LmrK>niE7 z^a1P-oWTtSNy)$uiHs`ePjCnbKw*bKARr%YYZ$PP=#$m5gw-GHC0hyxi`Ii~dq-=vHY>I$(t&gvr zZJD+_uHysi1|sCw5F9%@J2p1986V(Jz?-+`+Gu$E=DAFBo98{hwz#-LS7!Rw)4y6e{1mpkqVj;ek5QIr zo$y8x0+x0(kn!MkSy|catO<~2K~@6NtX+2m&ZN2aQ3osy#Bs;bv?Fkr)|K1<@yDAd zKeJrbvqkUn^J@`BOw%*QE8LFBAGS@Y*C`b|VJvNXOcMx=$>Embc0>uOnbTbXT@83^5)eoif0c$~k zr#O8+R6UXgKr}cwC^F}ZcF)Yr%-wx2MMJs?1Q@{C+ltS=CV=j1Nl8iYiUPT#vsK)q*a`Ra1oVV4+y2glo?y;t|66rSBY1+IZl zFFo5{RGlzn)1Uw#;>H4a7Q}-81&E2UCDCGjTj4A=waoEq>wa=jVqRT5jfy(n{wBf2 zg%`>mGc`4(at?O^!3Smxy@3D@QrFf*g|U~{L5|NH4ugPz=VYbnQI*7e&~Z*)Ufv}^ zh6g@1|I*3ckE3Pl{%6$qkznUbOG_Ykf#^gnDJ`9pdY#~S!-)B|YR;;hBZS5L3 z@bTkpdr(V|HGt$V($CMf8-F|lrZZ{w!T~$`4!SX0HPs>~aDD*+7H%yP8DMd#1;j#bg*X8A95w9LiEG6H_ zz`*InuTNX$FJ8R3y1MF#A`IE6X!tpjb3U%Ep=_KRrK_V;SX6X_R+w8%>@;3vE9e=k zH`f-RktNK-$q8HgQUv%U?H!BI@xEqL7nkiAaxPq4TmYZ~uRTE} z>FMr%&Sy?XLYa^2lLw7Z^5xP%&t`*yqpP8g zzJAB?TmYb`8(R(Y(y(zPqk%d&tOe|}>FenU+kbuuyllc{zSZv+pqn$`gUrQuLfng* z()@fY(2oI}B03=V`@k9`(X@7 zbl?5nBNcFzE$nRMJhlMr4xE+?yE~Ob+ofrl!Xf4eF621zwS{8lTgdfgaBy(BUPakF zf9CV&wyqU5Miq{em0)T%PfJ_7_`wUS4(oo04CGu@szfg(ATGaouiO8nzc4Vv%6W|I zH>c}V*E!)OX(u-sHA#`?Jy;!me-sRe0;nF=Kxe5qFL-^RjOHuwQh*PHJL*uL#dop| z(MvES1wuOv-s_`mYHA8O9WsYSCB{cbzYuXXQ&NI^VW4&bQ#%OKP)9&3I$a~6JiDLo z6HZ0B!0x1^Y6C_Hs}mO$ZG5Z>z~0`zBQ~-1gbFC@pyhnD{f!pbmzdZcFu9oo8TZBB63v25I|4#N zbHF6G;vxaNNdPPbWTj9PwC*0{6-98XtTs9%dehf2wWNaAyxm0N=3_Or;WRFNy1?Bo zMfxBB%tZhQO##~mF9&f>!DBQI){7|D@7V3eqyBwCKoUtxN*V?9-+89KufHEm>|Qpl z}+W{K)m+M!lYiKlRER|1gpm3YDA}OfQiw8r`dC-VxXmA$*`3lFp z5fdMO33{C3KHC={fMgEGwQXxo^R6%KoO8h5=EE@X^eT*WA3ogua62K6PD-4Ot$3iS z5bAN8M6*WfWFZ`aE?d%34b}iCcmrR^(UctotXikVrx#8Pc~?`jRSrrLa2lWtJVF5O z)7#&lo}NBQcik^@byrmME3D4!dK+LCy$~#Qta}@5BrY!Q%$S1Lgo=d30~lnQM>U-n zvF4>;X+kW&2bdb>9VC}*JK_O0cHO_vkW78GrnuyPv|F<=*LFUVBjISLfY>+$tn>=( zXsaAFm^nB&0HZCJl4fQu%+Fs0#sv)6!EL&>w$^6o`EZ_0$Y3f5NY%9=x12tLQlrIt ze+iS1kC(Sy8rwW@PhOTt8dwfMh=i!9FW_gY&$*+Y zJ$r0)UsF@l&+qimt{qfHTQDalr>%E@9mT|e(=my?u7gSgMSz*Q&gd2ynPBk62>`1q zv3@Lc^j1)ggNKNLGV^^rZ92E1mVrT|-!I!@q+id1hBz2GY3jn+h{|Wc2Gk7{Ry?Pa zL<1-R@g&M?AQOKa90WAo@k(_Z<_)TCkR04V+6)5tc@!@zD|^VcR-bnYLb=IzZ}EnD z9q_y9WxkG!FILseY-|mMflOjxPDgp|Q^?gX;LU(vwE%}$Utb3k0w5Sqf#Z4m?;@#( zWltnHO`Nb@aE>756c-oI_^`0C(M$Te0sIsb69f5EAQTWL^$I;ePW^%L6bkl$rvm!x zyYK{SJwGX^ zm0%QQt&5V^Btz8Q%ETmFuUrjWN=gQGu(Y5_FgGIXf3SKDoPKSxYW|I34a93}+8e+( zO+2fBfB-u?J1?&bIDZxvmP3t3K#huX{Q*5M3Aw(aiMIu%ZQZ9QL>yR(&lzDTzked3 zP*B=2dUX-t`98DhiGbY*89DjW$VgDDpPilxJ0?^AU!2lhDa58At?a#b`$Aeo2 zdO%qN1_zA1euCNrT^fL`Fms_e-4$Zq~_*cRfSw~v9c Date: Wed, 29 Nov 2023 11:25:53 -0500 Subject: [PATCH 14/18] refactor changepoint model --- src/gms/cp_gm_pb.jl | 119 +++++++++++++++++++++----------------------- test/gms/cp_gm.jl | 6 ++- 2 files changed, 62 insertions(+), 63 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 909d79e..57a5e54 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -25,7 +25,7 @@ struct NoEvent <: EventRelation end """ -Parameters for change point model +holds parameters for change point model """ struct CPParams <: GMParams # prior @@ -41,6 +41,9 @@ struct CPParams <: GMParams death_factor::Float64 end +""" +constructs parameter struct for change point model +""" function CPParams(client::Int64, objs::Vector{Int64}, mprior::MaterialPrior, pprior::PhysPrior, event_concepts::Vector{Type{<:EventRelation}}, @@ -142,7 +145,7 @@ update latents of a single element end """ -in case of collision: Gaussian drift update of mass and restitution# +in case of collision: Gaussian drift update of mass and restitution """ @gen function _collision_clause(pair_idx::Vector{Int64}, latents::Vector{BulletElemLatents}) latents[pair_idx[1]] = @trace(update_latents(latents[pair_idx[1]]), :new_latents_a) @@ -150,15 +153,15 @@ in case of collision: Gaussian drift update of mass and restitution# return latents end -@gen function _no_event_clause(pair_idx, latents::Vector{BulletElemLatents}) - return latents +function clause(::Type{Collision}) + _collision_clause end """ -Returns the generative function over latents for the event relation +in case of no event: no change """ -function clause(::Type{Collision}) - _collision_clause +@gen function _no_event_clause(pair_idx, latents::Vector{BulletElemLatents}) + return latents end function clause(::Type{NoEvent}) @@ -176,63 +179,60 @@ function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelati end end -# map mcmc kernel (link: https://www.gen.dev/docs/stable/ref/mcmc/#Composing-Stationary-Kernels) -# set of proposal functions - -# change randomly clause choices of mass -# revise for which event, modify the categorical -# e.g. another event type for the same pair or another event for one event type - -@gen function event_switch(clause, events, start_event_idx, pair, latents) - switch = Gen.Switch(map(clause, events)...) - return switch(start_event_idx, pair, latents) -end - """ -iterate over event concepts and evaluate predicates to active/deactive +map possible events to weight vector for birth decision using the predicates """ -@gen function event_kernel(state::CPState, params::CPParams) - obj_states = state.bullet_state.kinematics +function calculate_predicates(obj_states, event_concepts) object_pairs = collect(combinations(obj_states, 2)) - pair_idx = repeat(collect(combinations(1:length(obj_states), 2)), length(params.event_concepts)) + pair_idx = repeat(collect(combinations(1:length(obj_states), 2)), length(event_concepts)) pair_idx = [[0,0], pair_idx...] # for no event - # map possible events to weight vector for birth decision using the predicates - predicates = [predicate(event_type, a, b) for event_type in params.event_concepts for (a, b) in object_pairs] + predicates = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs] + events = vcat(NoEvent, [repeat([event_type], length(object_pairs)) for event_type in event_concepts]...) + return predicates, events, pair_idx +end - # birth of one or no random event - weights = copy(predicates) - active_events = copy(state.active_events) +function normalize_weights(weights, active_events) for idx in active_events # active events should not be born again weights[idx-1] = 0 end weights = [max(0, 1 - sum(weights)), weights...] # TODO: objects that are already involved in some events should not be involved in other event types as well - weights_normalized = weights ./ sum(weights) - - # Draw born event - events = vcat(NoEvent, [repeat([event_type], length(object_pairs)) for event_type in params.event_concepts]...) - - start_event_idx = @trace(categorical(weights_normalized), :start_event_idx) - if start_event_idx > 1 - push!(active_events, start_event_idx) - end - - - - updated_latents = @trace(event_switch(clause, events, start_event_idx, pair_idx[start_event_idx], state.bullet_state.latents), :event) - bullet_state = setproperties(state.bullet_state; latents = updated_latents) - - # death of one or no active event - weights = [(idx+1 in active_events && idx+1 != start_event_idx) ? max(1. - predicates[idx] * params.death_factor, 0.) : 0.0 for idx in 1:length(predicates)] # dying has a much lower chance of being born + return weights ./ sum(weights) +end +function calculate_death_weights(predicates, active_events, start_event_idx, death_factor) + can_die(idx) = idx+1 in active_events && idx+1 != start_event_idx + # dying has a much lower chance of being born + get_weight(idx) = can_die(idx) ? max(1. - predicates[idx] * death_factor, 0.) : 0.0 + weights = [get_weight(idx) for idx in 1:length(predicates)] weights = [max(0, 1-sum(weights)), weights...] # no event at index 1 - weights = weights ./ sum(weights) # normalize + return weights ./ sum(weights) +end - end_event_idx = @trace(categorical(weights), :end_event_idx) - if end_event_idx > 1 # nothing happens when no event dies - delete!(active_events, end_event_idx) - end +""" +create a Switch combinator to evaluate the corresponding clause for a started event to trace new latents +""" +@gen function event_switch(clause, events, start_event_idx, pair, bullet_state, active_events) + switch = Gen.Switch(map(clause, events)...) + updated_latents = @trace(switch(start_event_idx, pair, bullet_state.latents), :event) # trace latents + bullet_state = setproperties(bullet_state; latents = updated_latents) + start_event_idx > 1 && push!(active_events, start_event_idx) + return bullet_state, active_events +end + +""" +iterate over event concepts and evaluate predicates for newly activated events +""" +@gen function event_kernel(active_events, bullet_state, event_concepts, death_factor) + predicates, events, pair_idx = calculate_predicates(bullet_state.kinematics, event_concepts) + weights = normalize_weights(copy(predicates), active_events) + start_event_idx = @trace(categorical(weights), :start_event_idx) # up to one event is born + bullet_state, active_events = event_switch(clause, events, start_event_idx, pair_idx[start_event_idx], bullet_state, active_events) + + weights = calculate_death_weights(predicates, active_events, start_event_idx, death_factor) + end_event_idx = @trace(categorical(weights), :end_event_idx) # up to one active event dies + end_event_idx > 1 && delete!(active_events, end_event_idx) return active_events, bullet_state end @@ -241,25 +241,22 @@ end for one object, observe the noisy position in every dimension """ @gen function observe_position(k::RigidBodyState, noise::Float64) - pos = k.position - # add noise to position - obs = @trace(broadcasted_normal(pos, noise), :positions) - return obs + @trace(broadcasted_normal(k.position, noise), :positions) end """ -for one time step, run event and physics kernel +run event and physics kernel for one time step and observe noisy positions """ @gen function kernel(t::Int, prev_state::CPState, params::CPParams) - # event kernel - active_events, bullet_state = @trace(event_kernel(prev_state, params), :events) + active_events, bullet_state = @trace(event_kernel(prev_state.active_events, + prev_state.bullet_state, + params.event_concepts, + params.death_factor), :events) - # simulate physics for the next step based on the current state and observe positions bullet_state::BulletState = PhySMC.step(params.sim, bullet_state) - obs = @trace(Gen.Map(observe_position)(bullet_state.kinematics, Fill(params.obs_noise, params.n_objects)), :observe) + @trace(Gen.Map(observe_position)(bullet_state.kinematics, Fill(params.obs_noise, params.n_objects)), :observe) - next_state = CPState(bullet_state, active_events) - return next_state + return CPState(bullet_state, active_events) end """ diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 81f47be..7ee8637 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -31,7 +31,7 @@ function forward_test() trace, _ = Gen.generate(cp_model, (t, cp_params)); println("") - display(get_choices(trace)) + #display(get_choices(trace)) end function add_rectangle!(plt, xstart, xend, y; height=0.8, color=:blue) @@ -124,13 +124,15 @@ function update_test_2() # find first collision in the trace start_event_indices = [trace[:kernel=>i=>:events=>:start_event_idx] for i in 1:t] t1 = findfirst(x -> x == 2, start_event_indices) - +# TODO: validate existence of event # move first collision five steps earlier cm = choicemap(fixed_prior_cm) cm[:kernel => t1 => :events => :start_event_idx] = 1 cm[:kernel => t1 - 5 => :events => :start_event_idx] = 2 trace2, ls2, _... = Gen.update(trace, cm) trace3, delta_s, _... = Gen.regenerate(trace2, select(:kernel => t1 - 5 => :events => :event)) + # TODO: check ls2 + # print choices and @assert delta_s != -Inf @assert delta_s != NaN From 94f318f59f8ab22c368059e59d1e55e5893beea8 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 29 Nov 2023 12:08:27 -0500 Subject: [PATCH 15/18] test regenerating event --- src/gms/cp_gm_pb.jl | 2 +- test/gms/cp_gm.jl | 21 ++++++++++++++------- 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 57a5e54..3800e49 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -228,7 +228,7 @@ iterate over event concepts and evaluate predicates for newly activated events predicates, events, pair_idx = calculate_predicates(bullet_state.kinematics, event_concepts) weights = normalize_weights(copy(predicates), active_events) start_event_idx = @trace(categorical(weights), :start_event_idx) # up to one event is born - bullet_state, active_events = event_switch(clause, events, start_event_idx, pair_idx[start_event_idx], bullet_state, active_events) + bullet_state, active_events = @trace(event_switch(clause, events, start_event_idx, pair_idx[start_event_idx], bullet_state, active_events), :switch) weights = calculate_death_weights(predicates, active_events, start_event_idx, death_factor) end_event_idx = @trace(categorical(weights), :end_event_idx) # up to one active event dies diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 7ee8637..40f43d4 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -124,20 +124,27 @@ function update_test_2() # find first collision in the trace start_event_indices = [trace[:kernel=>i=>:events=>:start_event_idx] for i in 1:t] t1 = findfirst(x -> x == 2, start_event_indices) -# TODO: validate existence of event + # TODO: validate existence of event # move first collision five steps earlier cm = choicemap(fixed_prior_cm) cm[:kernel => t1 => :events => :start_event_idx] = 1 cm[:kernel => t1 - 5 => :events => :start_event_idx] = 2 trace2, ls2, _... = Gen.update(trace, cm) + @show ls2 + choices = get_choices(trace2) + display(get_submap(choices, :kernel => t1 => :events)) + display(get_submap(choices, :kernel => t1 -5 => :events)) + trace3, delta_s, _... = Gen.regenerate(trace2, select(:kernel => t1 - 5 => :events => :event)) - # TODO: check ls2 - # print choices and + @show delta_s + choices2 = get_choices(trace3) + display(get_submap(choices2, :kernel => t1 => :events)) + display(get_submap(choices2, :kernel => t1 -5 )) @assert delta_s != -Inf - @assert delta_s != NaN + @assert !isnan(delta_s) - return trace, trace2 + #return trace, trace2 end # redraw latents at same event start @@ -207,6 +214,6 @@ end #visualize_active_events() #constrained_test() #update_test() -#update_test_2() -update_test_3() +update_test_2() +#update_test_3() #switch_test_static() \ No newline at end of file From 8afaccc017508c40d84d41aa05a75acee4378be2 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 29 Nov 2023 15:32:04 -0500 Subject: [PATCH 16/18] WIP refactor switch --- src/gms/cp_gm_pb.jl | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 3800e49..0be3a9c 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -62,6 +62,9 @@ function CPParams(client::Int64, objs::Vector{Int64}, CPParams(mprior, pprior, event_concepts, sim, template, length(objs), obs_noise, death_factor) end +event_concepts = Type{<:EventRelation}[Collision] +switch = Gen.Switch(map(clause, event_concepts)...) + """ Current state of the change point model, simulation state and event state """ @@ -182,14 +185,15 @@ end """ map possible events to weight vector for birth decision using the predicates """ -function calculate_predicates(obj_states, event_concepts) +function calculate_predicates(obj_states) object_pairs = collect(combinations(obj_states, 2)) pair_idx = repeat(collect(combinations(1:length(obj_states), 2)), length(event_concepts)) pair_idx = [[0,0], pair_idx...] # for no event predicates = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs] - events = vcat(NoEvent, [repeat([event_type], length(object_pairs)) for event_type in event_concepts]...) - return predicates, events, pair_idx + event_ids = vcat(1, repeat(1:length(event_concepts), length(object_pairs))) + @show event_ids + return predicates, event_ids, pair_idx end function normalize_weights(weights, active_events) @@ -213,22 +217,21 @@ end """ create a Switch combinator to evaluate the corresponding clause for a started event to trace new latents """ -@gen function event_switch(clause, events, start_event_idx, pair, bullet_state, active_events) - switch = Gen.Switch(map(clause, events)...) - updated_latents = @trace(switch(start_event_idx, pair, bullet_state.latents), :event) # trace latents - bullet_state = setproperties(bullet_state; latents = updated_latents) - start_event_idx > 1 && push!(active_events, start_event_idx) - return bullet_state, active_events +@gen function event_switch(event_idx, start_event_idx, pair, bullet_state) + return end """ iterate over event concepts and evaluate predicates for newly activated events """ @gen function event_kernel(active_events, bullet_state, event_concepts, death_factor) - predicates, events, pair_idx = calculate_predicates(bullet_state.kinematics, event_concepts) + predicates, event_ids, pair_idx = calculate_predicates(bullet_state.kinematics, event_concepts) weights = normalize_weights(copy(predicates), active_events) start_event_idx = @trace(categorical(weights), :start_event_idx) # up to one event is born - bullet_state, active_events = @trace(event_switch(clause, events, start_event_idx, pair_idx[start_event_idx], bullet_state, active_events), :switch) + updated_latents = @trace(switch(event_idx, pair_idx[start_event_idx], bullet_state.latents), :event) + # TODO: use funcitonal collections + bullet_state = setproperties(bullet_state; latents = updated_latents) + start_event_idx > 1 && push!(active_events, start_event_idx) weights = calculate_death_weights(predicates, active_events, start_event_idx, death_factor) end_event_idx = @trace(categorical(weights), :end_event_idx) # up to one active event dies From cd870f2f7f684c572ef3fa157e767c1fc5e65683 Mon Sep 17 00:00:00 2001 From: Dominik Glandorf Date: Wed, 13 Dec 2023 13:08:41 -0500 Subject: [PATCH 17/18] refactor switch call --- src/gms/cp_gm_pb.jl | 29 ++++++++++++++--------------- test/gms/cp_gm.jl | 9 +++++++-- test/gms/plots/events.png | Bin 13144 -> 13463 bytes 3 files changed, 21 insertions(+), 17 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 0be3a9c..3969021 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -62,9 +62,6 @@ function CPParams(client::Int64, objs::Vector{Int64}, CPParams(mprior, pprior, event_concepts, sim, template, length(objs), obs_noise, death_factor) end -event_concepts = Type{<:EventRelation}[Collision] -switch = Gen.Switch(map(clause, event_concepts)...) - """ Current state of the change point model, simulation state and event state """ @@ -171,6 +168,10 @@ function clause(::Type{NoEvent}) _no_event_clause end + +event_concepts = Type{<:EventRelation}[NoEvent, Collision] +switch = Gen.Switch(map(clause, event_concepts)...) + """ """ @@ -188,11 +189,10 @@ map possible events to weight vector for birth decision using the predicates function calculate_predicates(obj_states) object_pairs = collect(combinations(obj_states, 2)) pair_idx = repeat(collect(combinations(1:length(obj_states), 2)), length(event_concepts)) - pair_idx = [[0,0], pair_idx...] # for no event + pair_idx = [[0,0], pair_idx...] # [0,0] for no event - predicates = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs] - event_ids = vcat(1, repeat(1:length(event_concepts), length(object_pairs))) - @show event_ids + predicates = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs if event_type != NoEvent] # NoEvent excluded and added in weights + event_ids = vcat(1, repeat(2:length(event_concepts), inner=length(object_pairs))) # 1 for NoEvent return predicates, event_ids, pair_idx end @@ -200,7 +200,7 @@ function normalize_weights(weights, active_events) for idx in active_events # active events should not be born again weights[idx-1] = 0 end - weights = [max(0, 1 - sum(weights)), weights...] + weights = [max(0, 1 - sum(weights)), weights...] # first element for NoEvent # TODO: objects that are already involved in some events should not be involved in other event types as well return weights ./ sum(weights) end @@ -217,18 +217,18 @@ end """ create a Switch combinator to evaluate the corresponding clause for a started event to trace new latents """ -@gen function event_switch(event_idx, start_event_idx, pair, bullet_state) - return +@gen function event_switch(event_idx, pair_idx, latents) + return (pair_idx, latents) end """ iterate over event concepts and evaluate predicates for newly activated events """ -@gen function event_kernel(active_events, bullet_state, event_concepts, death_factor) - predicates, event_ids, pair_idx = calculate_predicates(bullet_state.kinematics, event_concepts) +@gen function event_kernel(active_events, bullet_state, death_factor) + predicates, event_ids, pair_idx = calculate_predicates(bullet_state.kinematics) weights = normalize_weights(copy(predicates), active_events) start_event_idx = @trace(categorical(weights), :start_event_idx) # up to one event is born - updated_latents = @trace(switch(event_idx, pair_idx[start_event_idx], bullet_state.latents), :event) + updated_latents = @trace(switch(event_ids[start_event_idx], pair_idx[start_event_idx], bullet_state.latents), :event) # TODO: use funcitonal collections bullet_state = setproperties(bullet_state; latents = updated_latents) start_event_idx > 1 && push!(active_events, start_event_idx) @@ -252,8 +252,7 @@ run event and physics kernel for one time step and observe noisy positions """ @gen function kernel(t::Int, prev_state::CPState, params::CPParams) active_events, bullet_state = @trace(event_kernel(prev_state.active_events, - prev_state.bullet_state, - params.event_concepts, + prev_state.bullet_state, params.death_factor), :events) bullet_state::BulletState = PhySMC.step(params.sim, bullet_state) diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 40f43d4..297c13a 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -119,11 +119,16 @@ function update_test_2() cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) # generate initial trace - trace, _ = Gen.generate(cp_model, (t, cp_params)) + trace, ls = Gen.generate(cp_model, (t, cp_params)) # find first collision in the trace start_event_indices = [trace[:kernel=>i=>:events=>:start_event_idx] for i in 1:t] t1 = findfirst(x -> x == 2, start_event_indices) + @show ls + choices = get_choices(trace) + display(get_submap(choices, :kernel => t1 => :events)) + display(get_submap(choices, :kernel => t1 -5 => :events)) + # TODO: validate existence of event # move first collision five steps earlier cm = choicemap(fixed_prior_cm) @@ -139,7 +144,7 @@ function update_test_2() @show delta_s choices2 = get_choices(trace3) display(get_submap(choices2, :kernel => t1 => :events)) - display(get_submap(choices2, :kernel => t1 -5 )) + display(get_submap(choices2, :kernel => t1 -5 => :events )) @assert delta_s != -Inf @assert !isnan(delta_s) diff --git a/test/gms/plots/events.png b/test/gms/plots/events.png index f071e247da5431950ba16bad0bbea8030fe9c691..1924a8b8834e99a874bf3ea161ed7df79dad56a0 100644 GIT binary patch literal 13463 zcmaL8bzIcl_B}pyNlB@6cdC?xz<@Li(k-cgfJh7=Ee#?e9Xf(^N;fJEB8_zS(4F6d z&%O8a-1z;*KLlSh^FC*vz1LcM?PI9A>Jxn2d$OvqW@3B$Azg+5s z_JO}J&6J)%A=fv5J~m{>LLhVy1*nuJJQX?R?o4XgD0$^}7~?Fu#--BCgpQs6-j@qT zj)87T972i86F?b)Eg(GWv64McT|C_lYgt+m4iM$}C@lK9thg*Rf(Y#yPFX-1^IP)3a%$Vwi}CiSK?;Kw^f3gbe+lV`5;C zF(0M$Tki$D4+*&yDJ2hu2DSZp`%?*QlQ+Aq zdwWi;-Ro&g3f@g7NS<H|i)uc1Y`~t7xR5=a?#-ZfE!fs-()z{yXiOB(&Eo+)ZN-u3Mcq2Z= zB_y=}_>nUexO{Ty4I?eo-pHc0!ch@|3MkwT!j~~jjd^HY`-@Ds6qUDFr)-nCdJ4Wm zt(*#c4Wt6h;J|E78zqqK!9b)e}s~Tk!gxA<d7n z-XQy17XBezrtJhh6VwQ9LD|CNFKqL<+3yo`wC_XbBBbQ+F({Pis%w)YdQTw`1a{SO zpVqhg{6rXpSlvYbks4xh^7qWDK9}8&j;0?^l|_HUmqMI-dh9nE9)ewwLWlb_`5Nu9 zoYOgyjjJB+a_V%kBSR#VP(wWywl&57NsXw=$?}av#&5gZZ}HCU5TA*u3+CW$=~sO7eOI=Q9x$^_f&R&Z`pcNE?!$}`sAwx4lhOEBEK0$Td5YT6wkVJ%QD%oQ zW)BVB!tdnToQg^!UjA%s@QC*ZUFtLMg~2cIe_z(xhn{GGGX@14x5uQ7!n#4uxhj;I znH{~-EAAif?Cw5a(ahmql8<4+hVVq(iwvDps0#H5xpuGYCFNEVt+$UktwT)J6GlEE zDWb&@oSwq5#9AA^yPumhuS#_g}4~IX#JVoxt0a^AO$I`l4>oMJ&rN zpuuFqE_e#dYlUKTkuER&tty@GAgKXb0ZIKRfWlg^NUJAItg2-`n>v$pmv8Z+oC z3lkkZeK|gHATK8US^(207h9+-yWZedi=*ke>xD~Hj3pVpLVVM8!_|In7=y>r4?@Ye z^{fa!5b!A?*|-EwJw8X9dr??R9*-0vdn|=Qa#)bml=KV?xb%NV>vn7fA#@aeens+# z1Bg8o+Dx2lp(iVU`&aY)cLbZ;G&b=tfPOC`^q#bE1Ov0JanoJhsOg5?>=wm%t+%5a z$g%!K$Ds^3aw52z8sP|V-|Ib7Pom$LF0-xv)E!g{>jmW0&08m4rrmFseC9vyMRXSe zj#(x~T?1j!(!BYmijvu$xtoK5LvhO{96BUV#U=AE2KkF{cJI)y1$AZ~9N-Ae%2_a( zRS0=~LoQR{Ac&q3#YkYI9&OgZ-}Me?6c}+dk{YU^!d(riV==lv(ZJ=-|}kh;-j5;B)1T3$;$q8hGW(w#C|%) zHE@R}@9KS3mwSx3l3>9_HLuAMx!^dxevf=74eN_r*GLkQRcAh91-zH%2g&vy0sALI zXBsI0IQgiRn|lK)a>6!o<+fLEdNO==*MEcxd+B$cr4gWE{c04KbJi#;V7Tyx5h;Ij z$R<;hZ*Hf3%|Q=S2h>2zgmZlo{OzG9DujTxE$ln4j<$5+8i_XNv`2Wp8W9GDTjcFE z2*nyRF3sOFOG-WyO`tDN5t`84B9KYw@B1?irkQEcyQ=sR)%8*G6HYyC`rCAxCZVBz zQ$=Qfu}kglejf$z%Cgy>#mSp8+_oh&vXDsGeA1*scC{ELsmHxsRD2&{YRaVHf!Nf8u=BnUYO3f7ZP5Ot9rEpobEt?I#2Wnjli;wH0E5&^=6=4o(=G!#z1i!}Yq^A|ra+0JaN|-7msJn_{4(`WN|sPyLumcmp^$ zDCTz%Kzeb{`Jy^$A)A9_F;fZw0mto!^y6j=p@`=KIR5I_Plm59WhMz=!eCE9&RY=} zc=rxF9@_HRB~JEMU~FOT?UnGU+2_(#4~fWxhw-`CbcoX9zMg!jSnGTxW@bKYFD-8U zGiM3~cIVmu(0Bu!RR$nJY+Q~*Y-@5-EvABZMU5cq9XxaDK#EAo$0q6yd=I^*$qn>Mv)5N*Qapm*WQK$2ry$X@`ZT~m$zY(A} zs)-1F3<+Rt2bD&*Fyx(jSwXEB%>%xnxRCXBKtiu-N-UKEVj0qXXe=!i?YG8XHj^K0 zhANDLJS0Y}+RjXqk_|6<2@+sss}ViK4|xB;V%JRUfCzyA=VqmK2V-+OQVjqPBCkup z`rUQ^li{s^4#-UwAh5$RT3hMI4Yg*kIY@N^pNB_MkYT_y>u?y%XHD<^r#=RVE7M~I zZwQGZ>|eIyX~syw-wjLurb*00(PIU4kpLzp{R??*@_4QnLyXvKOj^s3q5@QRfr)X_8ZX9lSK`G02gB$op!>9!QM0ddF~ecf8iNB1FDu*NNBXq4>UOV z>>XsG#M~~Ow5QQur_|=TZA3-QsF27nW8cy!iJI~&pW+KaoPpviiHzcp%?~_}0`sGz z6D6owwp!?MQU8$HvvtkD*k+Dr!X+-8#xq|m^a6)G9FmpqNYaFc+B}pY{RzU#iqEYD zw?rKzVZ!n5Rnwfc*Fen96_NjEGyyd~tID79^MDdZN81QVk$wTn=|GA|9rkK>eeU~7 zx5?tWw6#*K4%yM7J2Ns?Ze1j@Fe3D)kK$C zr%sN6dw1|Qe+OWAlG-RKL>fi@vXR}M257?BZuB`*ClyB2eCQc!VfJselBr`@CFom> z7Oe~@tTrwVy(3-y*&(X2ZefW5*$YK#783o+o2(EAHOkZ4{@GRsfLxfN%ROPgA5e=- zj-hy`3-2a6Y=62TPX}#s<2U?2sEGr)LC6Lgy~b)tK71$(yM+Nsm5AINR6NKmQKGKg z@wodI8aDS2gd>i={hyMOEpG~iI1q!zir&a(l%P<*lltNKCYAczx9HB#Z!}RflL8TB zNYdtDIYnjyAEpbImVDuTuts8}nfY@BZDsQauCAY_ik9%2A7f!5Iwd@ZjVy&8*Wi}4 zMl74zwRddKa9okDsa|#%Zfb6I596@nK|oBaph+k3YCb{TcEdfG||7_wX`e?gQGwsX!90TdI%fL z+lRJ%*&-5cc;HroDZN&&-Q8=#@QEMZwk9QtX0VC8ce$3G@kYW#Uq4x~7xdMEP7=64^t7EBL3^aF@vFvk*v9KFrvguiaR5m4yy-*^U&z>aY;qu}lwazwAJgpu|-uH8m%WKwsMK zo%6%_vXQweyF~=4vuR{EJ4`NzRXu`8hcqJw2NFlnA9AI4RGa#0T_Yj#&Hk z+|0~UT&mE$;j{Xrgao%i6`>l`{BmO84?HfUqY|{!eevRjv2nWo)J+t03+W=l!pigu zTU%RgTmbixnpM0bwx6i%-8gibWqT1zM<<}s=}%5O8_9Hzz;Nr~qHS`{M_j-r#bZ zIvp$4w#{w{9Nlp@O;Cp1g0_6_;OBtB8ybS!WHy zRg=ZO73Eq7L$t;Y4V{VEd~M0DVipKF_*MqrHLLKXV_?w0k6IuBJ+q#DN-gD9wG6=4z3H&m^zK zy1G9UWYplFpo|&@1Usk$V%j5@=bvBbG$Le) z(5oi!o7bGZcUrK(n%s9=DV{!kI&2?-ec#$_kVUt!B`;&PP2O}x^}^hzL=YC-9^f|J zGFD1aBud4Xv~Hpy6VrpQVzGDnD+{meYrqI|bc^n75)wjM399Ot0~Ui>{yBOoXK7%$GDEdszGSc~d- z?_YUwE7~rWg`S%`Ou1=*ySm=oKmq=Ad_RZP+ImDSm>c_t6)fSg-jVbz;>c-IV8IO|(;W3#IKm zhJ|^128A(!4(f>6h4YAVrI^`$j8}P&KbEk4 z`TTJV5pA(62{b6vPZ9#fKiRymsR^_Nr?%<-GRtV+B!Qjm@mAum?9M*?phmCZE#Hp! zTvN3!`L74#xecoMcwb0C(HR(k^2ER^F_VqmV1-nn07`8@s(9P;*-Zpt(WnFdTPPi~ zb5YA{d0IeB9Q{hStft^DI3e9%ZS(*W)=>UpiqCTFRGC>uo&LA9X=efQ=X~UL>}Ps8 z0r)JM{zg^1ZnxdeRgU$utI40JFg7yU!0x)%qcrj5AQrmYy&%;&v$8LO@vNFJ2I`4o=GpOEC ze~RId#K90?`tLl??cFp^y^3$b&yWOv9Q1&R8=sy1(|Pxn7#1lbJqj}IMwp*o#?U=j ze-rHQOE7-tP%-l<&}#3!FCAx_Tl2qK&F%LERM)J>_0O)kV9b4MY7-NHb7~Ra1)DFo z?Ix!ZCH-+G3-Nuu;<#!EjIgYGls}1yg&k+-(V@V3@s42{5b~yGiV58YitYdGSNR!W zE6KBfz7*NjCH>+B(MEPUswy2bM#Iy9>JOX?$LE1`=r%8)u`go*%3IP79wVj7^0G z24A@e1&(XHV+&)6rNgF>$ZKntiE%66CJ9*n1Z74M9yEeE)zIM%6g45t> zzWf+?dvT!QNmBkY6Q3cJLQFB5!#ZyFd|E6dk~om|Mz&SHOngK^Ok5T=)^LYladCEP z3b!IFQ-8ZRkAsx#l?M8s{yL%F7%JPDCNXzfs2RBI2EJk^u%C*MK?UctIU}f570|x? z%V)T;LSpEe7yQ1{=gVheIWjV!zLO6e$4$=T;DIWWa)B(%LP6U(7<39Z{6DMaAd{_l zbfWZg=Nuna5m@4g^3=PKZ{#;2-M6YgN}gOJAr}gQ{QeEMO^|-c-o8-!Gif1%PwJoU z*WxiFgY?z`j@zt{km3~3o;f*eCzn8OTip8JZcR)pY87Q^X;~Pw7Db4GkvP9$jl{lc zvf^P#PV|X3BXYbFzozDJW@&$}lXYoiy9#5>7l;)~rmDDk;Q3CxOcpA>@~m@i@kROc zGKB5TO`;mZQ78q`cf+HH!#_QC6`Z!Cv)5u$h!3**70w_a@IIDgJSgNFz`647?cjE_ zC|2D=G+R}8L(b5}I*E&;6-Kl_=nVd)LH`<1EC3Yb2}Xm#q^P()BZpa;K3??;0Nd*~ zK0{t8@PI2!QTFyXnm8qCX??h=t{&Dr76iH4=i;Qq{&QtNMcKbo^I9UMW!!)#K)wm& zh>`1~@rwNGqExCkXY7-Es_|1))o2~Pt>@z@5-5qhuS{G8w;F^j(ue{@Zb~I=X#Pr> zVF0Ace~(X$oJ3*TUtCGn6;!n9&=7pzV4pLUa&tSn{QFq*<+Hyxqm(;D5vPInMC3;P zvdmSHI6O)_S!L8MTo;rjO~1+5I@f1NeN8RU&5%GGPx!Q{v=OXFJY-{^J zPuggGnp*M-zPmoCwh*h9Y`FM_84?#an_%i)Q-A4N5Bi!6Ev=j6W$-jHJ+#2Jk&_j%EAN>=nI|5O? zoxlbKX6WBf3qTksH(V?h5<5ZNoY&DfDQ5U^muzUy@-F z0M((Vw>nnZd00G!X(ts$17AF2krW!c@Xi=yFEH|)v$5myuyuHc1%sOc9i@`wuac51 z#*WxcuHnNStH8^kx;4UL!Pf_@PR?g9*m3cSEppizYPo>IOlS|-+oKyPMuk)x3XlG9 zr=W+~w&(P^?rx$$Gr9h!sawv;0J0r`g!yX}o-wdlD8!Q z_Sl|B;zC3><^89TC?l=F7lXeFlJt9t0`K3fbtWn=k8Wgi4`A4s-5;L@_o>+ChM_Dh(ItRl#{t&KzdYaVWCsaWnk7~4*W*o)C`fx6W@^TU#tIM#g}dQImv zrmVHDt+ch>+=qvwCubXY3BmZW2Kc^gVVpX+`PxQT(f+m+i{VM&8 z%KEja(-=qdoBb1nQD+%$G~^S16bX9f7tFf5e`X$2H_~j$sBUNV^Uk>?nR7ZO2Qfu( zIYoVsp<#a$o>R>NE)WRmN zW2H~s7!Z3QMMcsViV88Boc}U*9x>b)I|Eq@?J(`2JtsX`4%BO%S>?dn;_L+@ZY?pl zH8sK9&EFr0<6>pGv!$1&6zs(8&GJ4Yan!s zhU++UBJCXn>sJxA{aWUqX)&C8C(EIDyb|k&ivt zk=QwptQg{pc6~le2!Wwp?NMw_>RlTlLP%H_N4A@QT$*m+fBhaz44+Ssi-=J0)J?I3 z%jVGV+kT~IlqqZs`IfpTrK_dI>~yrfO$98a^704w?>|M#J+NNt`e5GvuJE`>!s}PI z+GXF_G9{STb;#?bymK~ZdmkHcZGtw)x=th;PKP}_L)_i{>WC}*_BW(!u|d_=(S(Bz zxPoyy7fgG0xIxj~T5H3}!7*OtFyo7Y#&gNH6i||rvpSUflv>nTy}qo7TfcnuWOw1# zH_>|IZQ84y;2{5<+ixVGa$RW-M8XdvCn7bt(=)$}94MqjvL>MTeh(&;Jf3`g>)XtE z)Af1NN`}wXa*Bg?9r4SXM@pQne(uJkCsceJ!nDuC289kmBQL&p|4kmV`$YWOK9)Qj zdiv_st1z9B3KO5p;Uy>JRGkh|u3R>T2}TO6f$_cY)kW@(Z-{$0MJIParENkid^>Jop|*2=3P~P|qC;{=Lup z7p+=(uVp2Wa#gu)PXk{hK0e+_cXVs2ZnQw}^}%X?W+qcgueR6u?m~Oh<i2h(|xW^Gcz;6w>f}d0?Wv7wvD&( zgn^%bxBm$(FwEeyepo|DzJm+z>(&0xbuJr%=4}xYUS|i}Gc!n6r5H^OjhN_Ye!H;} z@Eyg#e#Ro8Qm-3SPrPrf3WdHMT%!_oe%?D_;=c3&7*OIzW9DFw-iO0FXhP)V`(Io97l_cz)ECgK|z5E z&oQcV9WG9nDln?61a5t{T1JukXSOipc*XhIMlrJj2M^D?sHm9=+XO}YyoFz%WLRR^ zcf#kKYEo5Kz!5Vkcack?!-e)vEJ@A*zBdGUipEob00HoRdPARO2xD>fmuC~5Vq zQ_aA@06CW~oq~pxlwYqE-ed*HjIpt?vN8dM7}AM{2`MR_R~M&q--FTQm6fRgo0nO3 z$8i)rdg$iK!piC?Mgv(|TA~!Puiq%DgFxlvrad+al@fTgzAAkX+wZv_%c63V$c>h3 z16iPOT}9I%auIi;Gcqy~60(g0ICQHW=K=^NRu4Tu~G4Z}P*tcklYuc^qvfB_@JNi5#1-?($-1EHImgt>4=GE?WZW^MRwt1b~*$`NRy3 z`yw$UTQv<;y@6yVk+&a!hTG;CKk&Shh0?{`%S0&(wQa-6xg(Ccs3hLNCn{f7O}S2h z5*-*EY&u)#gs`fm+w~)44{M}vo!|{3WWzD=$Zr9bn5=#UmX;>&=_)6O=^M+c1|~NG zJBh$g1n&OGl2h(@u8mEjoBBdR=1(Q(y+wu2oE25QTIaP;5?1>v}G zCh|c_uSAAg?!2g$v9WhMwOa(G&=&B7TOD{|hxr!G2{~nDd%Nwl*DI3Ode(98*}@*G z;aisKf;Cqg7Ik%XrKVDax9%N|ny96Vofg$?on3lz?BtFJPud@+I5hq)HfbVelqAXI=%FEG0;*W_)Hh%z(rgaiJ-@IJdEc@8T{16*FIaN3XZBB}xYSg;X8QIzw{02# zA2U}+RnNW$Vc5S|UDB`{eSW+(g>rE)s0K0Q83pytMJt~it~s5|^-07BiS-K`HJ(Xa?68@)VXFk;^uwC#vy95wL~j}u=z;w-M;>txMP1_Qm; zFtQjXjH`Ld!GpEImEIH~?~9Y>b3tZi%aH;-NaNXhemkwt1(*j;*SxlNb#=A3x7XB+ z8}H3q7}hcA{rU4o@u2DeJR~I}dkDw_OaPBd+i5)4bF+tD^8NsM3aVBFL?RIQc}7J zgZ*Un$;k;A#*QP8z##zk(DISL2l5A5HA659b|DG&1O_wvAda{0@XZ!hP5aO1t8p3E zpGaPxs_E;G0|>n{{ZyD2S-0Kz=*@}ap;q>wp18j)&?1~VMb4sWfCEaRDG2d4?-fp4 z!nCvwLABFLdS3)l+TVF-6Ia`ii>R(XS%{M4VQ0U~XW9(XEYHCNAP)H8FTjq$G#(wk zv?(COCHNPy`%6GY~3zXI{vS2=dac) z#`bh?31q%{9Ueoibm|j%=evwg-c43Hn8+MoozGkYsCLTNF1)(j_6aV_7u7D%o!dV} z&30o&u)sDmC!lldH-bn4wha<1B{p_Hq3I&)jI<1pe4>D*qMcoV2R>@^qwRW-UBL3F zMZJ)dA3~23+moLiQgP9a$-L_q-;KY!`8r!6km zFJfwHqt{6z(Tq<7odpK+^6~)Q0ku4AxRLh-U?o~eR8-TZ8k961v_|qp;BA;A#+W65 z7P1FxT#@4-Q~4z%Bm@OrdlLD9-~rpc={U^H%_ENCgE2f?!4wky8+AEIt4?V1#V{yR3B$?qL61X38^0ZOCU4-Mr;2!bEL zD9xb8$x>0V3*_R}#fl`J|IUCC@1?g4FA#tUCf=t&cv@OnVUq7fhJ-+N*7J*1Qv`dL zoO+Y)V?8!Br1HJExIjU}WMXFaI!9K4w&N0%Rffkpd)36{Xwy|iZ66j5+!Snw*8Aj{ zvGFvhgme(Wv^1~Z#hK7KIR=DWitsB22w>A`_>{|NB#pSJqlg;4;K;R8(XnHz()$Qi92Ng{@Au9Vn(1 zK!q;rLjWi0Oq#s?iUaWVDr~fFroupX0nlsa2Ddae4#72xVvobdx~aVAs_r5_S(lBjALeX(`L4Lxws~QLTr{buY=5Nr;T6Hu=`m`68=r-JJz<>d#*AHdQKydJW%=bn`f+o&iji@lnM z#`@?pq=AdBxos8#tv`N50MG#EV&AiXbUn8%ZV>c5+5`fdTd&OGX!AD)21a*0cm8op zX=&-Gu_tnGshby8V#OAQ&;=b`YU*TMhA$fxusfoY?dx~jK$sHZ;vQVFv&80t>n6PQ zJKxX<#dg5=W7>#_G?cGXtd*}#MnQ3moEWvQx@l~HawH=oJCry9#Si4Jf~snz?#-19c@>q* zo$rKbZYmiP;+&itQ9c*u;M7@_6Codg)Z7>?7E0+2z{C#=4F#V83Y>$3BTvr}C@(ef zLrrC6We|%VtNj>erMwkCqPtjqTrxxGfArng--i2b%@@aQskY0v!S$S*dn47=)u7kv z=;#Qd(T2fxK|i9Ysd=Lj^YV*ql$QW0ri!|B%{c;ruC+3$W8(D?@F0*aU$U}5OM&D6 z3#4q4BxDPqP4!JCq!JC`le;Pp)WL2R> I(x(3Z54xnN!vFvP literal 13144 zcmaKTbzGF)*7hLM-3U?x0@9!~g22!Mf^;|1A>G|bH;4j)(%oGKBGTQUba#Gx^f~7| z?|IJe8~%_P@44@_*Iw&d*Sgjw{JDY@7CIR^1Oma5k%lQlAa^t%5QJFNJK&RR_3#1k z53-S*6by0;|9RU|kN|;DLu6p0DsGwEv+hoWrj5vZ1;exJPUikWxCyw_cVX09ajGr| zrbI}h6hsoZqE2^LH{4V+1z$h9k2KQo+PHL@5%V=Ai3!7JdTLy;Sa}l?RJyzPl_y%d zr!i*nVKUE#MpZ}8-^^3->yFOawX6?kG3XwvdLwHli^5=Zg}5c;)KKUUeqd}^SWquD z9TfVG6E!?c`X#0q6pByw|NA)K3%B)XNbddn_Yo;rE|n~vEqWzQKWJs*A5Xj3vQ<{b zCgf@&avyrMG+!CkRC{%;k={0stxDoN;(xJ z^G1D^m6M}|LaTSgT$}aE=w!^Q8css@jU_L*UD;YMvD=O*ZR_2!osh>$d$M^@ED5lF zI)$)ltz>b_Z3>6oZKtNYK5%=kM%;dZJ+FvTDsZPhDWpQH!&?`&v0!lTzud4-O0%x6 zL2Att6oOpbifB*Gf+>H@n>ubA8lJ5HRo%;#yg0tO8zfjbRO6Fiz_traem)Je{(U;TlLaTpN-Y4hidSr_Os@ zZH_bNL#|MuFOlQF-`?u&W}Yo(^YD{3&#Y5Iil}tn9ZK3f1JjC^v!c}dm0YkwLtK2k zxF4iJe+)z_>h-2b?!BJITSv!%J%1h%9XyUvG+G8V!&wOVa<@qHU2UQ<$6NPTAvAGh zq(Ov_ciO1V&yynSW|r(aR)*Uhs;$HQc`vkO9V%HKJg}I(8%GMRx^Hr?FZrR)p&6uo zeR#rNYw_z11{SxOx##8Scq0xJVULt1j`AKq5<#5SzYO`~sXQs76r~3_PnQkup<^Q> zdK`BqG0TocHJ}r|wRn-&&8Vq{ha8w0-#N?>HUf_qkMq?{sfaPxocB}rFvoqEwAJt` zqZriNG(Zu1k9qoZsT&B zVTuiv6};(up{UAMYaJeN=fUQ$1_Wvms?bsMptGXO$R;pEf-`OrKi9<|Bf5uf9;m&)A;HgzYqIM2KRiA0760Da%f~i_Hd3%?;=ouZlToH;sgNqo_x-ag z2!6(&%9j#lwTlJorq6CJhekeXor=;ykZIPXWzej75)(gE^T+HKpkT+quPGvW|Gq~4 zsVGc#k@LTi2%dtGeMq@V?`?uV`;j0cMqONn#{QDS{TL!Djo`3&=l<)jJwN;L?wwB$ zkO`6cUKvN95|ec0u6;;GX!5gJ8%9yYQLZe83Q;N+D1^M=#x#$9E$LLf)N>G3Bw4e9o-l@23vT6`Z#+*0F z==V-KnYBSgL6;A61WL3$)O4%vEgYoJwX_;~ZPniv(_QD}-|Dnxvn4hrvJ?dzIF<-f zTf&Z9t%JhCq$U22wZ=htLV*Z+cB>ppb=ql9{@jl>(^L7m7NEQ~U zkGg@yhi*a&=F%P>$EPk^Vg87t|H`VCy*stGSPLYyRphjOz|abj=NWtA(b@84Bv}G9 z%yRu^VPua$Xr+u}m{TiS+h6)LD(zq|sNL~5SHpc67HA6G{(njwaJxPpUb&Vb`b37x zK$YT%SUlmI;6Q16czV5;i%U_z)@^7wJ2UqE6y-FFV*)~dX1K3RK(^uB43o-dyS8hQ zBcs2UEacmqeV~Uf9=8Yli+lEYc`Cl@YQ(}9pXWK>_xq!@-&~)GcxMV%*ilh>A5Gab z>GEW(wtJJ)=b#HoN_Tu6;7<(;^ZLKw1Yq14xIN;vq(jA+JR|EvnN94lpl(qM!pD^E zcS+dQ*j|&Z+I5^T_#3~RCeEF@)Rx7on=5oifr%O6D%;1JdhPjla(hA`3U9in642Fc zQ{__jx(VR;Pq{<(;D46<4lcs`n8;3Re$48L#rO2E6Fbbr5A7fH@o1dxXt_S6qXW(z z?+iZMs`@ICfyI*%f1I&$$22TLtPcji$g6pb0v9b@RQjUEf%aq@!{B-v0n*SgkwV^H zGLqHdr>-~Hhz8yNpNInr7d0KO?JzW}0lq zp-n;a{qt2On~EJ4;;!>UVKQDW00Jevxudm0Zt-V1wI+XhHuf_;y}rYLu^IClNQ zm6P*at8nG#6y$|pS7^mmdAupP^$26J=Bw@6*&VKo{`6XLI`kc~_ix`>qi_6pFZY*g z{|UpuDX+=iFBotv^pM1*6qxkuO+KafIBL&rNo1Keu(05B-k7qfn~u%IBoKKE^Zzt( z4+%xHSN}Gmoa2F2t-WEu?ol{_Y_d)HY9F(mV3wg>+jf?)Zj)_TPzR{xi!RJwv-cPv z?W*e3OpnybA%7QG3rH^f8IW8!X&Ch7|&M zu;(D%+x5*YFTe0C)9dMpJT7%Vq%>dvwnp9fxQstlV~$TKYU|r{j;0{*^yYyvf$#jN zzuY0)P3J(E_lf-B;ZtE{_K9?1LDZ{waGgLxf|2p5iVC;uyn_31o;>_jt5hQ2K)SAt;@lKU4< z(8SK#7~{qOZ%8cFlaQ2LU9zcrr(NekvR+3q5tU)rOCU~UT|y})USTM{YBXKDHu)k$ zAq>_FCqE9flwLpPH6xqN;K;C;Jg3WG-9!roZBYscNE4? zN1q>dRj0cRuq^W!oZ2?GN@XBh8a__X&20e{hlm5*6&bNv&xVKf(&y{m*Qv+PF(M`4 z(9+St^&eY)LRl$qV8bs!z;iA9!8z zFKG-!0N=fyRh=^0M#kwqS|7Ka-gc>)EVdAY<$#LydTTavKsd1L-Ox`h1XsjRY{PY$ zq*{$8D&JyYm@<}YbbRZcHZI=BsNWp+U!KZi=aQ1ek89-92yvnlru6ej;7FbIBR6|w zeZj{u`7ZD;$PwUjVw9C-A|YY#EyzZ;f5oeF_d6{Q@kbei8QEFqtq*^2G2W*{H}cyz z$QTs@q}~lLjRY_=kj4A+w*rX2I=5FsWSKDv0g- z9orVw->6JYM=@D}n847rtYq=3waBwchmlRqB62+Zb??nJD#ozcCTyh(T*U853?iHP z1-NKCF=;VgOS}2|Xi@p}y&xDFAnijbxN&pw$nZ7c!JDwVpMDd>;9||HpMf~Sd#s7z zm6BpLjfWIIG!I={QSJ0cjH zLI45~i4d^OzR9Zc>CI~vUP$}0zpn44aaX7PITuP2{zsbb=cae-246lViv~R;84nD$ zkVXV!y6E{Z*y4HWA@pWt#UwhN1`1WyUghl^E>IkJ^Q6gq;LY|_&2TLP6bcDTcgs3(4`L$1!9Arc3(Ay!BxZ z&=h3gSy$6(`-@i$%6+w@zuy&#gG0{qZ8770G@v1nHq_2IS}9@@5-MRYjUwAH&<&hK zQdeu0fHwB++o%8f#5lB$XQk zKBu)GQTHARxbDnU7)xQ1=k$JB?d%BTp1b_Y{Cz@0uEsSfk@{^D1M=40+PAKD-^lUX zOCCpkTE54FnaP>;uPGm5i77sUH(cNQ};F_sIldqudcts5F)%fK;Zhk7?VnUaTq<^=I%(NlIFd;^N`Or9T59k!0H7yg8<+si{!#>RT!& zjh60t9NEWMTsNVYQ`HxqR%s|xCSq=G)F^#?e+oeEA+P14)0oRqVi=(bnV?H?->S#u zb{z&f`cBK~Txx3S-kyC)Q;wwRmWI~>TJzk-ZgdzUigoLBPme}h>p|(Xoz&BTRhenT zt*7IqZfD5=0p1DK&7o5^Q!D#_G2r0KQh;!Ba_+XA>TQmd=zFZc&LA2B=(&^anw>ytNdCmcL;{22 z1u;-Y#Kgs8CPf`vaC&9Q&F=&~~rJ81^OT%@cStHCU9=X|DE{ z)Y0-&GpE{kVNMwOP;mDw}cSFYJt15bdR z`Ang-#4H*HHCo(Tpof;`&c`C>(W%?O@jU|^0br&f6C>k~Mik(He!w)5up%WE1OC!L ze6tZ#FnGvu@bN#2`uak`Zk9L8Dx9rWmyXua^?e@MP?eL6@Go-p4LNVwY z5Zlye<>{a{X-!z4jdFTGmC4~9(-VsI^`ywQONxWTq~iID$M$zvwT+$%0&DvyTvg^M zIqw|)vKm8+EB_;k++L)B1^@w~kS*D35zNf{OcWt7RvHm-x&I;}f0Uy@Q${{g%(Yo^ z^-@O7NBOvnPn{l5^)hPKc3c1*1rbWaiTdzmP|(xK%>^};;V-Ms?(0fWc@ETve4n{o zww3GicQ;K9X%q;VbQ%m>FBN{*MdQDOpZ{WKpq2*&cK8YIwD&~Hjy?NJ=wNsVp4A)k z?`HwNASMy&`t#qyMr{W+MQ2f&UK0hr(e)aw#S%=O9EuNi9jF_$3{s2fw)2X-UuV(P zPxLiN9bX+gc+e^te*!jm!8Hl;S;PW`O??XRfd*w~P$7#3CuH@=9->y7TZ)4^qb@!acJVho+q4G`KYUH@<= zoHO_hYPhc^BF9`gHCWzr)MTL1Ymbb@O?`Y3`aHgRk_}H4=n^$8XHc$3%p_f16i+pD z-?+im6017xy50oEi|H?Zy6rqiJJIfeSFEEPYRu^9g>CWmU53&MC%o(7AK1?LX=|{^ zLkzo^NaSQ^oYBNIC5zNhi(VLwCG|&7`xUb^dq4kDd8~Y#8C$4RT7G-2zndFVIE@t- z3zzWd&`eEbiJcH2Gq#EAb|+{>H3ygkP;^Gn?frc-fb1NDr~HG16c_KEiD)E-J@heQ z95HAzI=6q$G{cw-fvT1z5d1=wE%bV$h>=~r(&u)&-3*vr5CSAGJlndS0FbMfD3 zGoHym(a7!hYbY&~6x+mXsUOjQ3KrO)^HZ&)nRhntkERjt(=dzfi`4vj_;fhj$Q)8P z1Qy?{SCT=EkLcQr`#0SZ9=@MYW!hMXHA*-#5Y^_<_zWwkWqk`2`Lddqc;aNwlxR)2 zrh*Ywd-Zy|4Tk7^ty)hh9UYlSc&(pWd~D+`+NdL#2&jJi?Ou+_B(j{0=;_Y5tzg#o z?^)4*nf!2h(!9cxH%bx`Yw9M1e%H`TXlc1@Hre!O1F2|qK3z3!mJ|W3B>EckbEGn2 z77^*=@rLuEmCO9X$s9u7|> z)!bwY3u>o5$a|pKYaA0yL_Bupr<8lA{`lWA2#P&37Ol+CC#q1dgD*8bsia+ zn8H_-P@#ATG(5b_&cVH zwr`F$BIfi)2AWE@NxtqLEdr9t5;Qh2=uI-i5rYfvvE^m7S|| zu!W8l;wXfKjGDJ_g@tI*JO*fwM$(F>3aRZ@QTRyr|7ght*@UH<$Ph+$DW8VUq-*mE zQb$L(Cyr3=-(vf|7aDYQZ&>p0HhHWJE}v24^oIYbGj2<0It@7{{kEb0R6tL;YGg6fcmj7$TB70JHja;>D|w6#u()?o$( zP#Fz{D33NkMwGEmZQH^-T`72{tNddRRFfPN(fL-@6pROKaA{04?RTFF0{tHnB;jLB z+Og+ji6rr)@_hDZ^aZq^ z@Dp+$6+?WO0}jiJ zeuE!C&uj$m>n%E}g$L~+K_??WPX^t~ZeTJ^c zsf5weFRT<<&hfbtgb24X*qFb^7tqQmY3V4?u!4uf!%%B@3JN9}EY0h>tj8MBawH=} zr3D;+Wt~<0M97p&$tiyJrT=&q$y$8iP3jZrFO|;b6$>3GMcLk87VT5}&6a(?Kd9Cq z6b_9KV*-MSl$cm-yz;efWfT&I0S99P|08?yaw`(O=@OI4#z$nFu53BE@T#8DZ>CxE zNCt>SK!_tSRiydAm_M3^*6e~AM+rp4{~|%~;%LoNBKSO=rOOB|j4nC=p>e#@RHLcZ zasRs@QR)0P?j9o0RCXI5>NM<7dRHHuycKQJpo2OCquPWoP=E7JD+4b3HQZzh^aZtz zoi!z7nmFWotHl2~*o9A*-sclKH5QQsw18(`mAMos zj&v`gdY`W{`<{sxk43LV4zJHtTMSdU&b)J`ZLBi>^xEh8LPA2qVm7Q*r@?V`Z7>%R z1-tm#@GHt<60Yj>*?B0WljP8s7ea8EhiuTB(WkBfgSpvO0VcGqW?{A)L+Vj1>&!6K= zdjajMq?4hQjQX%kfPwBIkTsyig6v!ShFxMip;4cUr-l^6_i}KJ=`9Kw1cNq{cIA5;LlR(5cx znBKM>D^3Gue_tOin69?SbtY*{$j-K_(So0*u!ZTHXx^Ov%1%f~sIi5)e$N*1?ZV;( zpa^bCQZQ^$H5r0Up>^4nhyy0HIn`^Lt+&5t2?8}{X?fW@vXaGaw$b_UXFpE{ zE05X_3~Ne$M-!k-6)0xg&NVmmn;n*1trujgmFj4!sEjl^TLpHM4NC-!oNQ0mSx-8i zAFQsdtVF24ZEp5j{hBm6If;gbMpUMntFES&oRpN5ot+&U``r0@uEle_OrKpw*~rMK z=j+hWkX`Hf&vt)ApfnizD2(LTRkBct_y`FKUWSm{qGDkkmpC=(OxN06?lkS9qobEz zZf(uA2-4CbK3(l&Wd5}rBLb{D|CaVh&}ogsxEB+N86~UX)vx@#ycySd@49&(&nmNi zVB;_CHk;{MpgCC#=O1rRyRQ9sM;RkEEnHLSbL~+zDP-6g6v6$XKw*32<8W?_@Dz}7 z6h*F8y5_16Ix*P0KYf}9w{x9!9tWz(FYR;SiQe8`;fuB0j~_qc2%o*Jr6?TOJf3qN z*VS!wN)vRo%@+1X#w1;E5-`WLH=7dJLC zl8ho%a=y;Kz512Rsx5pnVO*kB<8eHy?#+C7cnB6R=(07@wA(7k$+`b4+jlI-l$(%{ z@Ug{E_v!p$jL3~fh46z1@Acb!aGo4JFHp$3INjCP)7#O$1(wq>Fbrh!S4l(?M}qls zi3$tq6G}y(HK~{BvGzIo`1qKZtPJK#B{5Mt0om$-O#JSNam*7Qiy^xCZ5r_q7R^dn zvT10S!*VolsYZoZ6>-fAEv=~%?K;gWGu69`{)i~)k1ZTmzY=LCs~2m;#m6h}Z=e}1 z_eKNxu+44XXf!FF4#;v|=c_4|EF{2co2etP-DKP!q6rujf`Wq1=18$aY3b?Bz&a~g z+V*?N#X~UnJz6SRM9xy=b>n%ye$@>WW~z?Ccr2?>cOV|#D03xo>4Bb|^? zb9wpmSNOE6DUg*QSCKserC{!~y<;JPugemEe4 z5}S}fKtiJKP=bpv&*9+UAeX{M!D-YD(EjpdyUBjBb7VvX0`9}ARlWBl9|Ud+yS@wG z{dAYYw8Cr(754<4Q+$9mpL6sCO@T{(4P?Us!_udLuxacAos4RXBw4g%%%@zHeQH%L zEdh|&VB!)(&WnBP9?qoq06{+qa(`B+SG#@GRO^l|L(IACk4hM1oNw`Rvrz67@Z*Sudci_ zlf5&(^!06y7Nr39D!IF#*E=i=xNOy%e8HbjKRa_P8}|gjT3^ox0l_+sb!~anda+Ri zP!m)_jHfOxE?`H76-8QJ4(5P6M)bw?ReG%ty)O)$hsCC*l2TFGaiu{pDFluVR{L7( zKCSkpwCI)XvTNz-&B?=cnSw&C21iry4uG3E2PwULdbR8B+U)l@J2l=LKIIc;sD#(EC(EZuUoOgSTg&L?mC|U)O#y-1mN>bP#Cd z2op`N?2yo|3V^)x&EG9(e>w^5wqEXl;m7-e0`=xw5bRCPo3(~r$dF1F5QkwW__L0E ztYs}HQy{%Vaj51>b(?tW z%S(1SRYYV=KF1YlTU$22=H_OQvFB#SZ95OhxJ<%ATA1OPtZd-LjKYtm}GG|6=Gltlj-n2MQrcw7LJ zcQ0q?>FL=t$|3Cl#Ba`)V+OKCZo!b?^S5s^ZN6=zmZOcDHqBsv)6&xPy^oXZ9Du26 zMb0aqqz^2=O5)A+B>_G@n_lw-z;e&s7G5467m(;+ zjaMG}mo%8NJ{|O*}vR1&SU(07U=?p!~`uvs6r4vuf9Z=^N24T{PvB zi`RlXkiB+9oJ>B4r8k~C=f6a5`Ry0l!Hoza288$f!7t!>0M|f;wwWogAKRI>YkT5; zVA5ms?Afy()KWk~%Czg)=;`6$PwalcF<&|up-b%2bPb9*d`!n$H^XK2jnu|-7$z#( zAw}h_qN4H`JOZpwy;#4+;}~#Cty=5zAL2M@XLc`SWWIou0V_{PO?3k-R7J;n@#_3A z+v^~Fc)fYEw25p#`>YT?=m8iW(L-+6Jb0?E7oP_*j7q@i1HiP?IZwNucO77&3f%N8 z9*l1D!O&1NBVTp@3q$r<2=V2!EP+v@issv!OVyaDleNJv*Xr-f`z+eETL7H^%ehQ^ z_`=I`iw}zcbh7d@a$aDF1bQNir-N92+*;0)B- zVEWGWQY5?C_sq#k(?0grGqL{C(kc)OF`t&3ng#JjBaKLA=`D-tpQS!`csB=Vq5w1V zrti%!UnVvhauO0=tFht_9~igE?)tS|Z`%PLj476byud4#;{hV0ad$|5yq`|BBqt{) zHZ~Rvi6gkxOLZ)U^1|(s!#5KF0ZdHn?(J3cv{P+G6X@*fDpD9005P1-CTUqW9T^dk z4$6mfX;ERRUJLqL5F~Bq0ryg=Sn4JUu77dV`Qm7tWqNJn*Xgw1< z`+TLT9HTlhIePA^Nrj3cU~7N869@0{@$rp}jY(KEo)2$%V22v_M%@DpS6EmW zBnS=7qW|5xf9xY}(ytP{<^GWENHC=;U-X)9{Mf?6FRTk;JoS z-!6_f!Q;p&;IAiMXZC~iHxCXLOO>F$x~*OfOO6w4DMjVwdvQ{fX^+gm=SW0?5f;DQ zT=Tt0W+rB400G{D69ZH}XN!^St$Z(k0EU5#00O@o`~n8Tl>j+GWj@{<7rDJUL=IeB zP0_zO4!Joh5CPqWIryHTmR7nQbpFSW=a4c`4?wct$NmHW3IXtfyH9OxcL2cxlpiqt zO>hF->vYx?v_MhA>sOcaw|axWz{v(akP)XKhXKk0{&+S~Vf^UPBi}>jhf2n@1WH*? z4nEUDKwWIC1T+zl-jv-FP!7=yIs%Z0PUirKK!EeSA{7cA0p2#(?B3dQx!Wef$yrff zU%#|JAbe>JIK`*s^W!lo(2xM^up}o37z)-`SenLbrw4M+*4DNs9TfMp{-2$lVw%Y! zms^h@0GA5u$Lq{`B8b|~p}5>!SK~3#?yY-=7W1vwpbjvbNn&AP*$n6R1B&kzNzP-v zHq&spGuy;WN%fn-9RPw+ke8>UUJK3&0K)qQAT-D#uwlP01vC1&8ve z1PucN9M%L^mZY}WEp>+z#bk3B1`NLZ_We7hh>us7k5vve(2|w)M$_n7iRDa8Oh&Jx z-@SX+YV536#w#M?3jzkrj9Yis&(F`>usc)qJY?eY;111loePhO$^`oXdaVqwVjG+k zHeDt*wi{=F5bWwKy`Mg*&+ugQS5#C00#kNsCr7vf#E5jzG=U*&JpdSZZ_y|;XI2G@ zp%n%Rn-GX-K>A47bW&4N>Od}7Qr#yaaz5>OXcdLVNF5Ez69ASqt1@ALC1`RokB14wsaXJZ=oH)v^o!e*Rpy|Mf;M##K@1-k%u(gQR_e6*c#Dbxo!8rjM>~?3m0s}; zyi9Eevg_1y0I&gsB;m^lIM)IWLxAv}cO6eq6y^ub6T=&~96)S9H-Y%0)@B;?Cgjn* zpuGX-9_Z-kq9_DUfNNDT*fc9w(#+q09atMKh#7nd1P#+3cWT(uLXv6NhG?>B*|^nX zvZB(XtXsVltNxKC(077;wFD&4(Xq0}=;OK}07j6uufQpZBDE6R`o)k&6=C0-t6&Tg zC@n4M@tGZWK#XP!xgP*l4Bq2V30rfJf(o_utL!OZ9uVk0P+}r?IX^^gqy51EXfQwj zIWHA4CFR`57X==w3ZV4M$`-Ecfpvny?LkCmX>x``{YUDg%c-ue2273@XgiydVb$ Date: Mon, 18 Dec 2023 16:37:20 -0500 Subject: [PATCH 18/18] update tests --- src/gms/cp_gm_pb.jl | 35 +++++++++++++++++++++++---------- test/gms/cp_gm.jl | 48 ++++++++++++++++++++++++++++----------------- 2 files changed, 55 insertions(+), 28 deletions(-) diff --git a/src/gms/cp_gm_pb.jl b/src/gms/cp_gm_pb.jl index 3969021..38b03a3 100644 --- a/src/gms/cp_gm_pb.jl +++ b/src/gms/cp_gm_pb.jl @@ -125,10 +125,9 @@ function predicate(t::Type{Collision}, a::RigidBodyState, b::RigidBodyState) a_dim = a.aabb[2] - a.aabb[1] b_dim = b.aabb[2] - b.aabb[1] d = norm(Vector(a.position)-Vector(b.position))-norm((a_dim+b_dim)/2) # l2 distance - clamp(exp(-15d), 0., 1.) + clamp(exp(-15d), 1e-3, 1 - 1e-3) end -# gen functional collections """ update latents of a single element @@ -168,12 +167,11 @@ function clause(::Type{NoEvent}) _no_event_clause end - event_concepts = Type{<:EventRelation}[NoEvent, Collision] switch = Gen.Switch(map(clause, event_concepts)...) """ - +TODO: this function was intended to check if some event relations are impossible to be created a certain time step """ function valid_relations(state::CPState, event_concepts::Vector{Type{EventRelation}}) return event_concepts @@ -191,11 +189,15 @@ function calculate_predicates(obj_states) pair_idx = repeat(collect(combinations(1:length(obj_states), 2)), length(event_concepts)) pair_idx = [[0,0], pair_idx...] # [0,0] for no event + # break up to two lines predicates = [predicate(event_type, a, b) for event_type in event_concepts for (a, b) in object_pairs if event_type != NoEvent] # NoEvent excluded and added in weights event_ids = vcat(1, repeat(2:length(event_concepts), inner=length(object_pairs))) # 1 for NoEvent return predicates, event_ids, pair_idx end +""" +transform predicates for pairs of objects into a probability vector that adds to 1, including one weight for NoEvent at the first position +""" function normalize_weights(weights, active_events) for idx in active_events # active events should not be born again weights[idx-1] = 0 @@ -205,6 +207,9 @@ function normalize_weights(weights, active_events) return weights ./ sum(weights) end +""" +similar to normalize_weights but for death of events +""" function calculate_death_weights(predicates, active_events, start_event_idx, death_factor) can_die(idx) = idx+1 in active_events && idx+1 != start_event_idx # dying has a much lower chance of being born @@ -215,12 +220,22 @@ function calculate_death_weights(predicates, active_events, start_event_idx, dea end """ -create a Switch combinator to evaluate the corresponding clause for a started event to trace new latents +updates active events in a functional form +add=True -> add event to set of active events +add=False -> remove event from set of active events """ -@gen function event_switch(event_idx, pair_idx, latents) - return (pair_idx, latents) +function update_active_events(active_events::Set{Int64}, event_idx::Int64, add::Bool) + if event_idx == 1 + return active_events + end + if add + return union(active_events, Set([event_idx])) + else + return setdiff(active_events, Set([event_idx])) + end end + """ iterate over event concepts and evaluate predicates for newly activated events """ @@ -228,14 +243,14 @@ iterate over event concepts and evaluate predicates for newly activated events predicates, event_ids, pair_idx = calculate_predicates(bullet_state.kinematics) weights = normalize_weights(copy(predicates), active_events) start_event_idx = @trace(categorical(weights), :start_event_idx) # up to one event is born + updated_latents = @trace(switch(event_ids[start_event_idx], pair_idx[start_event_idx], bullet_state.latents), :event) - # TODO: use funcitonal collections bullet_state = setproperties(bullet_state; latents = updated_latents) - start_event_idx > 1 && push!(active_events, start_event_idx) + active_events = update_active_events(active_events, start_event_idx, true) weights = calculate_death_weights(predicates, active_events, start_event_idx, death_factor) end_event_idx = @trace(categorical(weights), :end_event_idx) # up to one active event dies - end_event_idx > 1 && delete!(active_events, end_event_idx) + active_events = update_active_events(active_events, end_event_idx, false) return active_events, bullet_state end diff --git a/test/gms/cp_gm.jl b/test/gms/cp_gm.jl index 297c13a..930ecaf 100644 --- a/test/gms/cp_gm.jl +++ b/test/gms/cp_gm.jl @@ -14,13 +14,13 @@ pprior = PhysPrior((3.0, 10.0), # mass (0.2, 1.0)) # restitution obs_noise = 0.05 -t = 120 +t = 80 fixed_prior_cm = Gen.choicemap() -fixed_prior_cm[:prior => :objects => 1 => :mass] = 2 -fixed_prior_cm[:prior => :objects => 2 => :mass] = 1 +fixed_prior_cm[:prior => :objects => 1 => :mass] = 2. +fixed_prior_cm[:prior => :objects => 2 => :mass] = 1. fixed_prior_cm[:prior => :objects => 1 => :friction] = 0.5 -fixed_prior_cm[:prior => :objects => 2 => :friction] = 1.2 +fixed_prior_cm[:prior => :objects => 2 => :friction] = 0.5 fixed_prior_cm[:prior => :objects => 1 => :restitution] = 0.2 fixed_prior_cm[:prior => :objects => 2 => :restitution] = 0.2 @@ -29,8 +29,8 @@ function forward_test() event_concepts = Type{<:EventRelation}[Collision] cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) - trace, _ = Gen.generate(cp_model, (t, cp_params)); - println("") + trace, weight = Gen.generate(cp_model, (t, cp_params)); + @show weight #display(get_choices(trace)) end @@ -83,10 +83,11 @@ function constrained_test() event_concepts = Type{<:EventRelation}[Collision] cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) - addr = 10 => :events => :start_event_idx - cm = Gen.choicemap(addr => 2) - trace, _ = Gen.generate(cp_model, (t, cp_params), cm) - display(get_choices(trace)) + #addr = 10 => :events => :start_event_idx + #cm = Gen.choicemap(addr => 2) + trace, weight = Gen.generate(cp_model, (t, cp_params), fixed_prior_cm) + @show weight + #display(get_choices(trace)) end # update priors @@ -119,7 +120,7 @@ function update_test_2() cp_params = CPParams(client, [a,b], mprior, pprior, event_concepts, obs_noise) # generate initial trace - trace, ls = Gen.generate(cp_model, (t, cp_params)) + trace, ls = Gen.generate(cp_model, (t, cp_params), fixed_prior_cm) # find first collision in the trace start_event_indices = [trace[:kernel=>i=>:events=>:start_event_idx] for i in 1:t] @@ -127,25 +128,36 @@ function update_test_2() @show ls choices = get_choices(trace) display(get_submap(choices, :kernel => t1 => :events)) - display(get_submap(choices, :kernel => t1 -5 => :events)) # TODO: validate existence of event # move first collision five steps earlier - cm = choicemap(fixed_prior_cm) + cm = choicemap() cm[:kernel => t1 => :events => :start_event_idx] = 1 cm[:kernel => t1 - 5 => :events => :start_event_idx] = 2 trace2, ls2, _... = Gen.update(trace, cm) - @show ls2 + #@show ls2 choices = get_choices(trace2) - display(get_submap(choices, :kernel => t1 => :events)) - display(get_submap(choices, :kernel => t1 -5 => :events)) + #display(get_submap(choices, :kernel => t1 => :events)) + #display(get_submap(choices, :kernel => t1 -5 => :events)) + + # the keys have to be enumerated, subsets do not work + trace3, delta_s, _... = Gen.regenerate(trace, select( + :kernel => t1 => :events => :event => :new_latents_a => :mass,)) + #:kernel => t1 => :events => :event => :new_latents_a => :restitution)) - trace3, delta_s, _... = Gen.regenerate(trace2, select(:kernel => t1 - 5 => :events => :event)) @show delta_s choices2 = get_choices(trace3) display(get_submap(choices2, :kernel => t1 => :events)) - display(get_submap(choices2, :kernel => t1 -5 => :events )) + + for i in 1:t + if project(trace3, select(:kernel => i)) == -Inf + @show i + display(get_submap(choices2, :kernel => i => :events)) + end + end + @show t1 + @show project(trace3, select(:kernel)) @assert delta_s != -Inf @assert !isnan(delta_s)