diff --git a/Project.toml b/Project.toml
index 538b118e9..b3e894556 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,7 +28,7 @@ MacroTools = "0.5"
 Parameters = "0.12"
 ReverseDiff = "1.4, 1.5"
 SpecialFunctions = "0.8, 0.9, 0.10, 1"
-julia = "1"
+julia = "1.3"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/docs/src/ref/extending.md b/docs/src/ref/extending.md
index 7f9dfd480..47f708768 100644
--- a/docs/src/ref/extending.md
+++ b/docs/src/ref/extending.md
@@ -97,6 +97,55 @@ apply_with_state
 update_with_state
 ```
 
+### Incremental computation for user-defined types
+
+Gen provides the generative functions [`Construct`](@ref) and
+[`GetField`](@ref) to support incremental computation for user-defined composite
+types, i.e., structs. They are useful when structs contain fields which are
+returned by other generative functions, or sampled from distributions. For
+example, we might define a struct representing 2D points:
+```julia
+struct Point2f0
+    x::Float64
+    y::Float64
+end
+```
+We might then use [`Construct`](@ref) and [`GetField`](@ref) within the
+following model of point motion:
+```julia
+@gen (static) function point_motion()
+    x ~ normal(0, 1)
+    y ~ normal(0, 1)
+    p ~ Construct(Point2f0)(x, y)
+    new_x ~ x_motion(p)
+    new_y ~ y_motion(p)
+end
+
+@gen (static) function x_motion(p::Point2f0)
+    x ~ GetField(Point2f0, :x)(p::Point2f0)
+    new_x = long_computation_for_x(x)
+    return new_x
+end
+
+@gen (static) function y_motion(p::Point2f0)
+    y ~ GetField(Point2f0, :y)(p::Point2f0)
+    new_y = long_computation_for_y(y)
+    return new_y
+end
+```
+Because [`Construct`](@ref) propagates changes to each field of `Point2f0`
+separately, writing the model in this way ensures that MCMC moves that adjust
+`x` only result in `x_motion` being re-run, while `y_motion` is not recomputed.
+
+Note that that `Construct` should be used with constructors `T(args...)`
+where the ``n``th argument corresponds to the ``n``th field of the type `T`.
+Default constructors meet this requirement. Other constructors are not
+guaranteed to propagate changes correctly.
+
+```@docs
+Construct
+GetField
+```
 
 ## [Custom distributions](@id custom_distributions)
 
diff --git a/docs/src/ref/modeling.md b/docs/src/ref/modeling.md
index 15be1990e..57e6e174b 100644
--- a/docs/src/ref/modeling.md
+++ b/docs/src/ref/modeling.md
@@ -461,7 +461,7 @@ The functions are also subject to the following restrictions:
 
 - Default argument values are not supported.
 
-- Julia closures are not allowed.
+- Constructing named or anonymous Julia functions (and closures) is not allowed.
 
 - List comprehensions with internal `@trace` calls are not allowed.
 
diff --git a/src/dsl/dsl.jl b/src/dsl/dsl.jl
index 37664742a..77e794663 100644
--- a/src/dsl/dsl.jl
+++ b/src/dsl/dsl.jl
@@ -147,7 +147,7 @@ function parse_gen_function(ast, annotations, __module__)
     return_type = get(def, :rtype, :Any)
     static = DSL_STATIC_ANNOTATION in annotations
     if static
-        make_static_gen_function(name, args, body, return_type, annotations)
+        make_static_gen_function(name, args, body, return_type, annotations, __module__)
     else
         args = map(a -> resolve_grad_arg(a, __module__), args)
         make_dynamic_gen_function(name, args, body, return_type, annotations)
diff --git a/src/dsl/static.jl b/src/dsl/static.jl
index cd00dee34..ba57c635e 100644
--- a/src/dsl/static.jl
+++ b/src/dsl/static.jl
@@ -58,7 +58,7 @@ gen_node_name(arg::Symbol) = gensym(arg)
 gen_node_name(arg::QuoteNode) = gensym(string(arg.value))
 
 "Parse @trace expression and add corresponding node to IR."
-function parse_trace_expr!(stmts, bindings, fn, args, addr)
+function parse_trace_expr!(stmts, bindings, fn, args, addr, __module__)
     expr_s = "@trace($fn($(join(args, ", "))), $addr)"
     name = gen_node_name(addr) # Each @trace node is named after its address
     node = gen_node_name(addr) # Generate a variable name for the StaticIRNode
@@ -91,7 +91,7 @@ function parse_trace_expr!(stmts, bindings, fn, args, addr)
         # Create Julia node for each argument to gen_fn_or_dist
         arg_name = gen_node_name(arg_expr)
         push!(inputs, parse_julia_expr!(stmts, bindings,
-            arg_name, arg_expr, QuoteNode(Any)))
+            arg_name, arg_expr, QuoteNode(Any), __module__))
     end
     # Add addr node (a GenerativeFunctionCallNode or RandomChoiceNode)
     push!(stmts, :($(esc(node)) = add_addr_node!(
@@ -101,14 +101,26 @@ function parse_trace_expr!(stmts, bindings, fn, args, addr)
     return name
 end
 
+function set_module_for_global_constants(expr, bindings, __module__)
+    expr = MacroTools.postwalk(expr) do e
+        if MacroTools.@capture(e, var_Symbol) && !haskey(bindings, var) && var != :end
+            :($__module__.$var)
+        else
+            e
+        end
+    end
+    return expr
+end
+
 "Parse a Julia expression and add a corresponding node to the IR."
 function parse_julia_expr!(stmts, bindings, name::Symbol, expr::Expr,
-                           typ::Union{Symbol,Expr,QuoteNode})
+                           typ::Union{Symbol,Expr,QuoteNode}, __module__)
     resolved = resolve_symbols(bindings, expr)
     inputs = collect(resolved)
-    input_vars = map((x) -> esc(x[1]), inputs)
+    input_vars = map((x) -> x[1], inputs)
     input_nodes = map((x) -> esc(x[2]), inputs)
-    fn = Expr(:function, Expr(:tuple, input_vars...), esc(expr))
+    expr = set_module_for_global_constants(expr, bindings, __module__)
+    fn = Expr(:function, Expr(:tuple, input_vars...), expr)
     node = gensym(name)
     push!(stmts, :($(esc(node)) = add_julia_node!(
         builder, $fn, inputs=[$(input_nodes...)], name=$(QuoteNode(name)),
@@ -117,17 +129,17 @@ function parse_julia_expr!(stmts, bindings, name::Symbol, expr::Expr,
 end
 
 function parse_julia_expr!(stmts, bindings, name::Symbol, var::Symbol,
-                           typ::Union{Symbol,Expr,QuoteNode})
+                           typ::Union{Symbol,Expr,QuoteNode}, __module__)
     if haskey(bindings, var)
         # Use the existing node instead of creating a new one
         return bindings[var]
     end
-    node = parse_julia_expr!(stmts, bindings, name, Expr(:block, var), typ)
+    node = parse_julia_expr!(stmts, bindings, name, :($__module__.$var), typ, __module__)
     return node
 end
 
 function parse_julia_expr!(stmts, bindings, name::Symbol, var::QuoteNode,
-                           typ::Union{Symbol,Expr,QuoteNode})
+                           typ::Union{Symbol,Expr,QuoteNode}, __module__)
     fn = Expr(:function, Expr(:tuple), var)
     node = gensym(name)
     push!(stmts, :($(esc(node)) = add_julia_node!(
@@ -137,7 +149,7 @@ function parse_julia_expr!(stmts, bindings, name::Symbol, var::QuoteNode,
 end
 
 function parse_julia_expr!(stmts, bindings, name::Symbol, value,
-                           typ::Union{Symbol,Expr,QuoteNode})
+                           typ::Union{Symbol,Expr,QuoteNode}, __module__)
     fn = Expr(:function, Expr(:tuple), QuoteNode(value))
     node = gensym(name)
     push!(stmts, :($(esc(node)) = add_julia_node!(
@@ -159,23 +171,23 @@ function parse_param_line!(stmts::Vector{Expr}, bindings, name::Symbol, typ)
 end
 
 "Parse assignments and add corresponding nodes for the right-hand-side."
-function parse_assignment_line!(stmts, bindings, lhs, rhs)
+function parse_assignment_line!(stmts, bindings, lhs, rhs, __module__)
     if isa(lhs, Expr) && lhs.head == :tuple
         # Recursively handle tuple assignments
         name, typ = gen_node_name(rhs), QuoteNode(Any)
-        node = parse_julia_expr!(stmts, bindings, name, rhs, typ)
+        node = parse_julia_expr!(stmts, bindings, name, rhs, typ, __module__)
         bindings[name] = node
         for (i, lhs_i) in enumerate(lhs.args)
             # Assign lhs[i] = rhs[i]
             rhs_i = :($name[$i])
-            parse_assignment_line!(stmts, bindings, lhs_i, rhs_i)
+            parse_assignment_line!(stmts, bindings, lhs_i, rhs_i, __module__)
         end
     else
         # Handle single variable assignment (base case)
         (name::Symbol, typ) = parse_typed_var(lhs)
         # Generate new node name if name is already bound
         node_name = haskey(bindings, name) ? gensym(name) : name
-        node = parse_julia_expr!(stmts, bindings, node_name, rhs, typ)
+        node = parse_julia_expr!(stmts, bindings, node_name, rhs, typ, __module__)
         # Old bindings are overwritten with new nodes
         bindings[name] = node
     end
@@ -184,32 +196,32 @@ function parse_assignment_line!(stmts, bindings, lhs, rhs)
 end
 
 "Parse a return line and add corresponding return node."
-function parse_return_line!(stmts, bindings, expr)
+function parse_return_line!(stmts, bindings, expr, __module__)
     name, typ = gensym("return"), QuoteNode(Any)
-    node = parse_julia_expr!(stmts, bindings, name, expr, typ)
+    node = parse_julia_expr!(stmts, bindings, name, expr, typ, __module__)
     bindings[name] = node
     push!(stmts, :(set_return_node!(builder, $(esc(node)))))
     return Expr(:return, expr)
 end
 
 "Parse and rewrite expression if it matches an @trace call."
-function parse_and_rewrite_trace!(stmts, bindings, expr)
+function parse_and_rewrite_trace!(stmts, bindings, expr, __module__)
     if MacroTools.@capture(expr, e_gentrace)
         # Parse "@trace(f(xs...), addr)" and return fresh variable
         call, addr = expr.args
         if addr == nothing static_dsl_syntax_error(expr, "Address required.") end
         fn, args = call.args[1], call.args[2:end]
-        parse_trace_expr!(stmts, bindings, fn, args, something(addr))
+        parse_trace_expr!(stmts, bindings, fn, args, something(addr), __module__)
     else
         expr # Return expression unmodified
     end
 end
 
 "Parse line (i.e. top-level expression) of a static Gen function body."
-function parse_static_dsl_line!(stmts, bindings, line)
+function parse_static_dsl_line!(stmts, bindings, line, __module__)
     # Walk each line bottom-up, parsing and rewriting :gentrace expressions
     rewritten = MacroTools.postwalk(
-        e -> parse_and_rewrite_trace!(stmts, bindings, e), line)
+        e -> parse_and_rewrite_trace!(stmts, bindings, e, __module__), line)
     # If line is a top-level @trace call, we are done
     if MacroTools.@capture(line, e_gentrace) return end
     # Match and parse any other top-level expressions
@@ -220,10 +232,10 @@ function parse_static_dsl_line!(stmts, bindings, line)
         parse_param_line!(stmts, bindings, name, typ)
     elseif MacroTools.@capture(line, lhs_ = rhs_)
         # Parse "lhs = rhs"
-        parse_assignment_line!(stmts, bindings, lhs, rhs)
+        parse_assignment_line!(stmts, bindings, lhs, rhs, __module__)
     elseif MacroTools.@capture(line, return expr_)
         # Parse "return expr"
-        parse_return_line!(stmts, bindings, expr)
+        parse_return_line!(stmts, bindings, expr, __module__)
     elseif line isa LineNumberNode
         # Skip line number nodes
     else
@@ -234,18 +246,18 @@ end
 
 "Parse static Gen function body line by line."
 function parse_static_dsl_function_body!(
-    stmts::Vector{Expr}, bindings::Dict{Symbol,Symbol}, expr)
+    stmts::Vector{Expr}, bindings::Dict{Symbol,Symbol}, expr, __module__)
     # TODO: Use line number nodes to improve error messages in generated code
     if !(isa(expr, Expr) && expr.head == :block)
         static_dsl_syntax_error(expr)
     end
     for line in expr.args
-        parse_static_dsl_line!(stmts, bindings, line)
+        parse_static_dsl_line!(stmts, bindings, line, __module__)
     end
 end
 
 "Generates the code that builds the IR of a static Gen function."
-function make_static_gen_function(name, args, body, return_type, annotations)
+function make_static_gen_function(name, args, body, return_type, annotations, __module__)
     # Construct the builder for the intermediate representation (IR)
     stmts = Expr[]
     push!(stmts, :(bindings = Dict{Symbol, StaticIRNode}()))
@@ -265,7 +277,7 @@ function make_static_gen_function(name, args, body, return_type, annotations)
         bindings[arg.name] = node
     end
     # Parse function body and add corresponding nodes to the IR
-    parse_static_dsl_function_body!(stmts, bindings, body)
+    parse_static_dsl_function_body!(stmts, bindings, body, __module__)
     push!(stmts, :(ir = build_ir(builder)))
     expr = gensym("gen_fn_defn")
     # Handle function annotations (caching Julia nodes by default)
diff --git a/src/dynamic/update.jl b/src/dynamic/update.jl
index bf4a42de2..3e3605f59 100644
--- a/src/dynamic/update.jl
+++ b/src/dynamic/update.jl
@@ -86,7 +86,7 @@ function traceat(state::GFUpdateState, gen_fn::GenerativeFunction{T,U},
     if has_previous
         prev_call = get_call(state.prev_trace, key)
         prev_subtrace = prev_call.subtrace
-        get_gen_fn(prev_subtrace) === gen_fn || gen_fn_changed_error(key)
+        get_gen_fn(prev_subtrace) == gen_fn || gen_fn_changed_error(key)
         (subtrace, weight, _, discard) = update(prev_subtrace,
             args, map((_) -> UnknownChange(), args), constraints)
     else
diff --git a/src/modeling_library/distributions/normal.jl b/src/modeling_library/distributions/normal.jl
index 53c90b947..c99bdcb7c 100644
--- a/src/modeling_library/distributions/normal.jl
+++ b/src/modeling_library/distributions/normal.jl
@@ -54,9 +54,8 @@ Float64
 const broadcasted_normal = BroadcastedNormal()
 
 function logpdf(::Normal, x::Real, mu::Real, std::Real)
-    var = std * std
-    diff = x - mu
-    -(diff * diff)/ (2.0 * var) - 0.5 * log(2.0 * pi * var)
+    z = (x - mu) / std
+    - (abs2(z) + log(2π))/2 - log(std)
 end
 
 function logpdf(::BroadcastedNormal,
@@ -65,17 +64,17 @@ function logpdf(::BroadcastedNormal,
                 std::Union{AbstractArray{<:Real}, Real})
     assert_has_shape(x, broadcast_shapes_or_crash(mu, std);
                      msg="Shape of `x` does not agree with the sample space")
+    z = (x .- mu) ./ std
     var = std .* std
     diff = x .- mu
-    sum(-(diff .* diff) ./ (2.0 * var) .- 0.5 * log.(2.0 * pi * var))
+    sum(- (abs2.(z) .+ log(2π)) / 2 .- log.(std))
 end
 
 function logpdf_grad(::Normal, x::Real, mu::Real, std::Real)
-    precision = 1. / (std * std)
-    diff = mu - x
-    deriv_x = diff * precision
+    z = (x - mu) / std
+    deriv_x = - z / std
     deriv_mu = -deriv_x
-    deriv_std = -1. / std + (diff * diff) / (std * std * std)
+    deriv_std = -1. / std + abs2(z) / std
     (deriv_x, deriv_mu, deriv_std)
 end
 
@@ -85,11 +84,10 @@ function logpdf_grad(::BroadcastedNormal,
                      std::Union{AbstractArray{<:Real}, Real})
     assert_has_shape(x, broadcast_shapes_or_crash(mu, std);
                      msg="Shape of `x` does not agree with the sample space")
-    precision = 1.0 ./ (std .* std)
-    diff = mu .- x
-    deriv_x = sum(diff .* precision)
-    deriv_mu = sum(-deriv_x)
-    deriv_std = sum(-1.0 ./ std .+ (diff .* diff) ./ (std .* std .* std))
+    z = (x .- mu) ./ std
+    deriv_x = sum(- z ./ std)
+    deriv_mu = -deriv_x
+    deriv_std = sum(-1. ./ std .+ abs2.(z) ./ std)
     (deriv_x, deriv_mu, deriv_std)
 end
 
diff --git a/src/modeling_library/modeling_library.jl b/src/modeling_library/modeling_library.jl
index 13d6e4880..34e621f62 100644
--- a/src/modeling_library/modeling_library.jl
+++ b/src/modeling_library/modeling_library.jl
@@ -78,7 +78,8 @@ include("recurse/recurse.jl")
 include("switch/switch.jl")
 
 #############################################################
-# abstractions for constructing custom generative functions #
+# custom deterministic generative functions #
 #############################################################
 
 include("custom_determ.jl")
+include("structs.jl")
diff --git a/src/modeling_library/structs.jl b/src/modeling_library/structs.jl
new file mode 100644
index 000000000..b298235da
--- /dev/null
+++ b/src/modeling_library/structs.jl
@@ -0,0 +1,103 @@
+export StructDiff, Construct, GetField
+
+"Represents differences in the fields of user-defined structs."
+struct StructDiff{T, D <: Tuple{Vararg{<:Diff}}} <: Diff
+    diffs::D
+    function StructDiff{T}(diffs::D) where {T, D}
+        return new{T, D}(diffs)
+    end
+end
+
+StructDiff{T}(diffs::Diff...) where {T} = StructDiff{T}(diffs)
+
+function get_field_diff(diff::StructDiff, fieldname::Symbol)
+    return _static_get_field_diff(diff, Val(fieldname))
+end
+
+@generated function _static_get_field_diff(diff::StructDiff{T}, ::Val{F}) where {T, F}
+    idx = findfirst(fieldnames(T) .== F)
+    return :(diff.diffs[$idx])
+end
+
+"""
+    gen_fn = Construct(T::Type)
+
+Constructs instances of a user-defined composite type `T` while supporting
+incremental computation. Changes to fields are propagated as long as `gen_fn`
+called with `args` invokes a type constructor `T(args...)` such that the
+``n``th argument corresponds to the ``n``th field of the type.
+"""
+struct Construct{T} <: CustomUpdateGF{T, T} end
+
+Construct(type::Type) = Construct{type}()
+
+function apply_with_state(::Construct{T}, args) where {T}
+    obj = T(args...)
+    return (obj, obj)
+end
+
+function update_with_state(::Construct{T}, obj, args,
+                           argdiffs::Tuple{Vararg{NoChange}}) where {T}
+    return (obj, obj, NoChange())
+end
+
+function update_with_state(::Construct{T}, obj, args,
+                           argdiffs::Tuple{Vararg{UnknownChange}}) where {T}
+    new_obj = T(args...)
+    return (new_obj, new_obj, UnknownChange())
+end
+
+@generated function update_with_state(::Construct{T}, obj, args,
+                                      argdiffs::Tuple) where {T}
+    atypes, ftypes = fieldtypes(args), fieldtypes(T)
+    constructed_by_field = length(atypes) == length(ftypes) &&
+        hasmethod(T, atypes) && all(map(zip(atypes, ftypes)) do (AT, FT)
+            return AT <: FT || hasmethod(convert, (Type{FT}, AT))
+        end)
+    if constructed_by_field
+        return quote
+            new_obj = T(args...)
+            return (new_obj, new_obj, StructDiff{T}(argdiffs))
+        end
+    else
+        return quote
+            new_obj = T(args...)
+            return (new_obj, new_obj, UnknownChange())
+        end
+    end
+end
+
+"""
+    gen_fn = GetField(T::Type, fieldname::Symbol)
+
+Returns the field value of a user-defined composite type `T` while supporting
+incremental computation. Changes to fields are propagated accordingly as long
+`gen_fn` is called on an instance of `T` created by [`Construct`](@ref)
+according to the documented requiments.
+"""
+struct GetField{T, F, FT} <: CustomUpdateGF{FT, Nothing}
+    function GetField{T, F}() where {T, F}
+        return new{T, F, fieldtype(T, F)}()
+    end
+end
+
+GetField(type::Type, fieldname::Symbol) = GetField{type, fieldname}()
+
+function apply_with_state(::GetField{T, F}, (obj,)::Tuple{T}) where {T, F}
+    return (getfield(obj, F), nothing)
+end
+
+function update_with_state(::GetField{T, F}, _, (obj,),
+                           (diff,)::Tuple{NoChange}) where {T, F}
+    return (nothing, getfield(obj, F), NoChange())
+end
+
+function update_with_state(::GetField{T, F}, _, (obj,),
+                           (diff,)::Tuple{StructDiff{T}}) where {T, F}
+    return (nothing, getfield(obj, F), get_field_diff(diff, F))
+end
+
+function update_with_state(::GetField{T, F}, _, (obj,),
+                           (diff,)::Tuple{UnknownChange()}) where {T, F}
+    return (nothing, getfield(obj, F), UnknownChange())
+end
diff --git a/test/dsl/static_dsl.jl b/test/dsl/static_dsl.jl
index e062d06f9..b602169a7 100644
--- a/test/dsl/static_dsl.jl
+++ b/test/dsl/static_dsl.jl
@@ -83,6 +83,28 @@ end
 @load_generated_functions
 end
 
+# for testing that macros can generate SML function definitions
+macro make_foo_sml()
+    return quote
+        @gen (static) function foo_in_macro()
+            a = 1
+            return a + 1
+        end
+    end
+end
+
+# for testing that Julia expressions inside SML functions can
+# depend on variables defined in external lexical scope
+module MyModuleC
+using Gen: @gen, @load_generated_functions
+external_var = 1
+@gen (static) function uses_external()
+    return external_var + 1
+end
+
+@load_generated_functions()
+end
+
 @testset "static DSL" begin
 
 function get_node_by_name(ir, name::Symbol)
@@ -603,4 +625,14 @@ ch = get_choices(tr)
 @test length(get_submaps_shallow(ch)) == 1
 end
 
+@testset "returning a SML function from macro" begin
+@make_foo_sml()
+Gen.load_generated_functions()
+@test get_retval(simulate(foo_in_macro, ())) == 2
+end # @testset
+
+@testset "global variables" begin
+@test MyModuleC.uses_external() == 2
+end
+
 end # @testset "static DSL"
diff --git a/test/modeling_library/distributions.jl b/test/modeling_library/distributions.jl
index 319668e8c..2c3cfb37b 100644
--- a/test/modeling_library/distributions.jl
+++ b/test/modeling_library/distributions.jl
@@ -104,6 +104,10 @@ end
     # random
     x = normal(0, 1)
 
+    # does not overflow
+    @test logpdf(normal, 1e13, 0, 1) == -5e25
+    @test logpdf_grad(normal, 1e13, 0, 1) == (-1e13, 1e13, 1e26)
+
     # logpdf_grad
     f = (x, mu, std) -> logpdf(normal, x, mu, std)
     args = (0.4, 0.2, 0.3)
diff --git a/test/modeling_library/modeling_library.jl b/test/modeling_library/modeling_library.jl
index e5e8d95ae..40edf391f 100644
--- a/test/modeling_library/modeling_library.jl
+++ b/test/modeling_library/modeling_library.jl
@@ -1,4 +1,5 @@
 include("custom_determ.jl")
+include("structs.jl")
 include("distributions.jl")
 include("choice_at.jl")
 include("call_at.jl")
diff --git a/test/modeling_library/structs.jl b/test/modeling_library/structs.jl
new file mode 100644
index 000000000..543946e9d
--- /dev/null
+++ b/test/modeling_library/structs.jl
@@ -0,0 +1,119 @@
+@testset "incremental computation for user-defined structs" begin
+
+    struct Point2D{T <: Real}
+        x::T
+        y::T
+    end
+
+    Point2D(x) = Point2D(x, x)
+
+    @testset "diff structs" begin
+        diff = StructDiff{Point2D{Float64}}(NoChange(), UnknownChange())
+        @test Gen.get_field_diff(diff, :x) == NoChange()
+        @test Gen.get_field_diff(diff, :y) == UnknownChange()
+    end
+
+    @testset "type construction" begin
+        gen_fn = Construct(Point2D)
+        trace = simulate(gen_fn, (1, 2))
+        @test get_retval(trace) == Point2D(1, 2)
+
+        # Expect NoChange() if no fields change
+        new_trace, _, retdiff, _ =
+            update(trace, (1, 2), (NoChange(), NoChange()), choicemap())
+        @test get_retval(new_trace) == Point2D(1, 2)
+        @test retdiff == NoChange()
+
+        # Expect UnknownChange() if all fields have unknown changes
+        new_trace, _, retdiff, _ =
+            update(trace, (4, 3), (UnknownChange(), UnknownChange()), choicemap())
+        @test get_retval(new_trace) == Point2D(4, 3)
+        @test retdiff == UnknownChange()
+
+        # Expect a StructDiff if fields change
+        new_trace, _, retdiff, _ =
+            update(trace, (4, 3), (IntDiff(3), UnknownChange()), choicemap())
+        @test get_retval(new_trace) == Point2D(4, 3)
+        @test retdiff == StructDiff{Point2D}(IntDiff(3), UnknownChange())
+
+        # Expect UnknownChange() if a non-default constructor is used
+        new_trace, _, retdiff, _ =
+            update(trace, (0,), (UnknownChange(),), choicemap())
+        @test get_retval(new_trace) == Point2D(0, 0)
+        @test retdiff == UnknownChange()
+    end
+
+    @testset "field access" begin
+        point = Point2D(1, 2)
+
+        get_x = GetField(Point2D, :x)
+        get_y = GetField(Point2D, :y)
+        trace_x = simulate(get_x, (point,))
+        trace_y = simulate(get_y, (point,))
+        @test get_retval(trace_x) == 1
+        @test get_retval(trace_y) == 2
+
+        # Expect NoChange() if no fields change
+        new_trace_x, _, retdiff_x, _ =
+            update(trace_x, (point,), (NoChange(),), choicemap())
+        new_trace_y, _, retdiff_y, _ =
+            update(trace_y, (point,), (NoChange(),), choicemap())
+        @test get_retval.((new_trace_x, new_trace_y)) == (1, 2)
+        @test (retdiff_x, retdiff_y) == (NoChange(), NoChange())
+
+        # Expect StructDiff changes to propagate
+        new_point = Point2D(0, 4)
+        diff = StructDiff{Point2D}(UnknownChange(), IntDiff(2))
+        new_trace_x, _, retdiff_x, _ =
+            update(trace_x, (new_point,), (diff,), choicemap())
+        new_trace_y, _, retdiff_y, _ =
+            update(trace_y, (new_point,), (diff,), choicemap())
+        @test get_retval.((new_trace_x, new_trace_y)) == (0, 4)
+        @test (retdiff_x, retdiff_y) == (UnknownChange(), IntDiff(2))
+
+        # Expect UnknownChange() if argdiff is UnknownChange()
+        new_trace_x, _, retdiff_x, _ =
+            update(trace_x, (new_point,), (UnknownChange(),), choicemap())
+        new_trace_y, _, retdiff_y, _ =
+            update(trace_y, (new_point,), (UnknownChange(),), choicemap())
+        @test get_retval.((new_trace_x, new_trace_y)) == (0, 4)
+        @test (retdiff_x, retdiff_y) == (UnknownChange(), UnknownChange())
+    end
+
+    @testset "test function" begin
+        @gen (static,diffs) function foo()
+            x ~ normal(0, 1)
+            y ~ normal(0, 1)
+            p ~ Construct(Point2D{Float64})(x, y)
+            a ~ GetField(Point2D{Float64}, :x)(p)
+            b ~ GetField(Point2D{Float64}, :y)(p)
+            q ~ Construct(Point2D{Float64})(a, b)
+            return p
+        end
+
+        Gen.load_generated_functions(@__MODULE__)
+
+        trace, _ = generate(foo, (), choicemap(:x => 1, :y => 2))
+        @test trace[:p] == Point2D{Float64}(1, 2)
+        @test trace[:a] == 1
+        @test trace[:b] == 2
+        @test get_retval(trace) == Point2D{Float64}(1, 2)
+
+        # Expect NoChange() if choices are not updated
+        new_trace, _, retdiff, _ = update(trace, (), (), choicemap())
+        @test get_retval(new_trace) == Point2D{Float64}(1, 2)
+        @test retdiff == NoChange()
+
+        # Expect changes to each field to propagate separately
+        new_trace, _, retdiff, _ = update(trace, (), (), choicemap(:y => 1))
+        @test get_retval(new_trace) == Point2D{Float64}(1, 1)
+        @test retdiff == StructDiff{Point2D{Float64}}(NoChange(), UnknownChange())
+
+        # Expect UnknownChange() if all fields have unknown changes
+        new_trace, _, retdiff, _ =
+            update(trace, (), (), choicemap(:x => 0, :y => 1))
+        @test get_retval(new_trace) == Point2D{Float64}(0, 1)
+        @test retdiff == UnknownChange()
+    end
+
+end
diff --git a/test/modeling_library/unfold.jl b/test/modeling_library/unfold.jl
index 3d8c1b54a..a77e2d256 100644
--- a/test/modeling_library/unfold.jl
+++ b/test/modeling_library/unfold.jl
@@ -1,15 +1,15 @@
-@testset "unfold combinator" begin
+std = 1.0
 
-    std = 1.0
+@gen (static) function kernel(t::Int, x_prev::Float64, (grad)(alpha::Float64), (grad)(beta::Float64))
+    x = @trace(normal(x_prev * alpha + beta, std), :x)
+    return x
+end
 
-    @gen (static) function kernel(t::Int, x_prev::Float64, (grad)(alpha::Float64), (grad)(beta::Float64))
-        x = @trace(normal(x_prev * alpha + beta, std), :x)
-        return x
-    end
+foo = Unfold(kernel)
 
-    Gen.load_generated_functions()
+@load_generated_functions()
 
-    foo = Unfold(kernel)
+@testset "unfold combinator" begin
 
     @testset "Julia call" begin
         @test length(foo(5, 0., 1.0, 1.0)) == 5