1
1
using ModelingToolkit, JumpProcesses, LinearAlgebra, NonlinearSolve, Optimization,
2
2
OptimizationOptimJL, OrdinaryDiffEq, RecursiveArrayTools, SciMLBase,
3
3
SteadyStateDiffEq, StochasticDiffEq, SymbolicIndexingInterface,
4
- DiffeqCallbacks , Test
4
+ DiffEqCallbacks , Test
5
5
using ModelingToolkit: t_nounits as t, D_nounits as D
6
6
7
7
# Sets rnd number.
536
536
end
537
537
SymbolicIndexingInterface. symbolic_container (s:: NumSymbolCache ) = s. sc
538
538
function SymbolicIndexingInterface. is_observed (s:: NumSymbolCache , x)
539
- return symbolic_type (x) != NotSymbolic () && ! is_variable (s, x) && ! is_parameter (s, x) && ! is_independent_variable (s, x)
539
+ return symbolic_type (x) != NotSymbolic () && ! is_variable (s, x) &&
540
+ ! is_parameter (s, x) && ! is_independent_variable (s, x)
540
541
end
541
542
function SymbolicIndexingInterface. observed (s:: NumSymbolCache , x)
542
543
res = ModelingToolkit. build_function (x,
584
585
end
585
586
end
586
587
end
587
- function SymbolicIndexingInterface. with_updated_parameter_timeseries_values (p:: Vector{Float64} , args... )
588
+ function SymbolicIndexingInterface. with_updated_parameter_timeseries_values (
589
+ :: NumSymbolCache , p:: Vector{Float64} , args... )
588
590
for (idx, buf) in args
589
591
if idx == 1
590
592
p[1 : 2 ] .= buf
604
606
dea2 = DiffEqArray (Vector{Float64}[], Float64[])
605
607
return ParameterTimeseriesCollection ((dea1, dea2), deepcopy (ps))
606
608
end
607
- function SciMLBase. get_saveable_values (p:: Vector{Float64} , tsidx)
609
+ function SciMLBase. get_saveable_values (:: NumSymbolCache , p:: Vector{Float64} , tsidx)
608
610
if tsidx == 1
609
611
return p[1 : 2 ]
610
612
else
@@ -614,38 +616,96 @@ end
614
616
615
617
@variables x (t) ud1 (t) ud2 (t) xd1 (t) xd2 (t)
616
618
@parameters kp
617
- sc = SymbolCache ([x], Dict (ud1 => 1 , xd1 => 2 , ud2 => 3 , xd2 => 4 , kp => 5 ), t; timeseries_parameters = Dict (ud1 => ParameterTimeseriesIndex (1 , 1 ), xd1 => ParameterTimeseriesIndex (1 , 2 ), ud2 => ParameterTimeseriesIndex (2 , 1 ), xd2 => ParameterTimeseriesIndex (2 , 2 )))
619
+ sc = SymbolCache ([x],
620
+ Dict (ud1 => 1 , xd1 => 2 , ud2 => 3 , xd2 => 4 , kp => 5 ),
621
+ t;
622
+ timeseries_parameters = Dict (
623
+ ud1 => ParameterTimeseriesIndex (1 , 1 ), xd1 => ParameterTimeseriesIndex (1 , 2 ),
624
+ ud2 => ParameterTimeseriesIndex (2 , 1 ), xd2 => ParameterTimeseriesIndex (2 , 2 )))
618
625
sys = NumSymbolCache (sc)
619
626
620
627
function f! (du, u, p, t)
621
628
du .= u .* t .+ p[5 ] * sum (u)
622
629
end
623
630
fn = ODEFunction (f!; sys = sys)
624
631
prob = ODEProblem (fn, [1.0 ], (0.0 , 1.0 ), [1.0 , 2.0 , 3.0 , 4.0 , 5.0 ])
625
- cb1 = PeriodicCallback (0.1 ; initial_affect = true , final_affect = true , save_positions = (false , false )) do integ
632
+ cb1 = PeriodicCallback (0.1 ; initial_affect = true , final_affect = true ,
633
+ save_positions = (false , false )) do integ
626
634
integ. p[1 : 2 ] .+ = exp (- integ. t)
627
635
SciMLBase. save_discretes! (integ, 1 )
628
636
end
629
637
function affect2! (integ)
630
638
integ. p[3 : 4 ] .+ = only (integ. u)
631
639
SciMLBase. save_discretes! (integ, 2 )
632
640
end
633
- cb2 = DiscreteCallback ((args... ) -> true , affect2!, save_positions = (false , false ), initialize = (c, u, t, integ) -> affect2! (integ))
641
+ cb2 = DiscreteCallback ((args... ) -> true , affect2!, save_positions = (false , false ),
642
+ initialize = (c, u, t, integ) -> affect2! (integ))
634
643
sol = solve (deepcopy (prob), Tsit5 (); callback = CallbackSet (cb1, cb2))
635
644
636
645
ud1val = getindex .(sol. discretes. collection[1 ]. u, 1 )
637
646
xd1val = getindex .(sol. discretes. collection[1 ]. u, 2 )
638
647
ud2val = getindex .(sol. discretes. collection[2 ]. u, 1 )
639
648
xd2val = getindex .(sol. discretes. collection[2 ]. u, 2 )
640
649
641
- for (sym, timeseries_index, val, buffer, isobs, check_inference) in [
642
- (ud1, 1 , ud1val, zeros (length (ud1val)), false , true )
643
- ([ud1, xd1], 1 , vcat .(ud1val, xd1val), map (_ -> zeros (2 ), ud1val), false , true )
644
- ((ud2, xd2), 2 , tuple .(ud2val, xd2val), map (_ -> zeros (2 ), ud2val), false , true )
645
- (ud2 + xd2, 2 , ud2val .+ xd2val, zeros (length (ud2val)), true , true )
646
- ([ud2 + xd2, ud2 * xd2], 2 , vcat .(ud2val .+ xd2val, ud2val .* xd2val), map (_ -> zeros (2 ), ud2val), true , true )
647
- ((ud1 + xd1, ud1 * xd1), 1 , tuple .(ud1val .+ xd1val, ud1val .* xd1val), map (_ -> zeros (2 ), ud1val), true , true )
648
- ]
650
+ for (sym, timeseries_index, val, buffer, isobs, check_inference) in [(ud1,
651
+ 1 ,
652
+ ud1val,
653
+ zeros (length (ud1val)),
654
+ false ,
655
+ true )
656
+ ([ud1, xd1],
657
+ 1 ,
658
+ vcat .(ud1val,
659
+ xd1val),
660
+ map (
661
+ _ -> zeros (2 ),
662
+ ud1val),
663
+ false ,
664
+ true )
665
+ ((ud2, xd2),
666
+ 2 ,
667
+ tuple .(ud2val,
668
+ xd2val),
669
+ map (
670
+ _ -> zeros (2 ),
671
+ ud2val),
672
+ false ,
673
+ true )
674
+ (ud2 + xd2,
675
+ 2 ,
676
+ ud2val .+
677
+ xd2val,
678
+ zeros (length (ud2val)),
679
+ true ,
680
+ true )
681
+ (
682
+ [ud2 + xd2,
683
+ ud2 * xd2],
684
+ 2 ,
685
+ vcat .(
686
+ ud2val .+
687
+ xd2val,
688
+ ud2val .*
689
+ xd2val),
690
+ map (
691
+ _ -> zeros (2 ),
692
+ ud2val),
693
+ true ,
694
+ true )
695
+ (
696
+ (ud1 + xd1,
697
+ ud1 * xd1),
698
+ 1 ,
699
+ tuple .(
700
+ ud1val .+
701
+ xd1val,
702
+ ud1val .*
703
+ xd1val),
704
+ map (
705
+ _ -> zeros (2 ),
706
+ ud1val),
707
+ true ,
708
+ true )]
649
709
getter = getp (sys, sym)
650
710
if check_inference
651
711
@inferred getter (sol)
678
738
end
679
739
end
680
740
681
- for subidx in [1 , CartesianIndex (2 ), :, rand (Bool, length (val)), rand (eachindex (val), 4 ), 2 : 5 ]
741
+ for subidx in [
742
+ 1 , CartesianIndex (2 ), :, rand (Bool, length (val)), rand (eachindex (val), 4 ), 2 : 5 ]
682
743
if check_inference
683
744
@inferred getter (sol, subidx)
684
745
if ! isa (val[subidx], Number)
706
767
(ud2, xd1, xd2),
707
768
ud1 + ud2,
708
769
[ud1 + ud2, ud1 * xd1],
709
- (ud1 + ud2, ud1 * xd1),
710
-
711
- ]
770
+ (ud1 + ud2, ud1 * xd1)]
712
771
getter = getp (sys, sym)
713
772
@test_throws Exception getter (sol)
714
773
@test_throws Exception getter ([], sol)
732
791
((kp, x), true , tuple .(kpval, xval), false ),
733
792
(2 ud2, true , 2 .* ud2val, true ),
734
793
([kp, 2 ud1], true , vcat .(kpval, 2 .* ud1val), false ),
735
- ((kp, 2 ud1), true , tuple .(kpval, 2 .* ud1val), false ),
794
+ ((kp, 2 ud1), true , tuple .(kpval, 2 .* ud1val), false )
736
795
]
737
796
getter = getu (sys, sym)
738
797
if check_inference
741
800
@test getter (sol) == val
742
801
reference = val_is_timeseries ? val : xval
743
802
for subidx in [
744
- 1 , CartesianIndex (2 ), : , rand (Bool, length (reference)),
745
- rand (eachindex (reference), 4 ), 2 : 6
803
+ 1 , CartesianIndex (2 ), :, rand (Bool, length (reference)),
804
+ rand (eachindex (reference), 4 ), 2 : 6
746
805
]
747
806
if check_inference
748
807
@inferred getter (sol, subidx)
767
826
((x, ud1), (_xval, _ud1val), true ),
768
827
(x + ud2, _xval + _ud2val, true ),
769
828
([2 x, 3 xd1], [2_ xval, 3_ xd1val], true ),
770
- ((2 x, 3 xd2), (2_ xval, 3_ xd2val), true ),
829
+ ((2 x, 3 xd2), (2_ xval, 3_ xd2val), true )
771
830
]
772
831
getter = getu (sys, sym)
773
832
@test_throws Exception getter (sol)
781
840
@test getter (integ) == val
782
841
end
783
842
784
- xinterp = sol (0.1 : 0.1 : 0.3 , idxs= x)
785
- xinterp2 = sol (sol. discretes. collection[2 ]. t[2 : 4 ], idxs= x)
843
+ xinterp = sol (0.1 : 0.1 : 0.3 , idxs = x)
844
+ xinterp2 = sol (sol. discretes. collection[2 ]. t[2 : 4 ], idxs = x)
786
845
ud1interp = ud1val[2 : 4 ]
787
846
ud2interp = ud2val[2 : 4 ]
788
847
@@ -796,13 +855,12 @@ end
796
855
(x, c2[2 ], xinterp2[1 ]),
797
856
(x, c2[2 : 4 ], xinterp2),
798
857
([x, ud2], c2[2 ], [xinterp2[1 ], ud2interp[1 ]]),
799
- ([x, ud2], c2[2 : 4 ], vcat .(xinterp2, ud2interp)),
858
+ ([x, ud2], c2[2 : 4 ], vcat .(xinterp2, ud2interp))
800
859
]
801
- res = sol (t, idxs= sym)
860
+ res = sol (t, idxs = sym)
802
861
if res isa DiffEqArray
803
862
res = res. u
804
863
end
805
864
@test res == val
806
- end
865
+ end
807
866
end
808
-
0 commit comments