diff --git a/examples/multi-objective/mdp/power/power-timed.nm b/examples/multi-objective/mdp/power/power-timed.nm deleted file mode 100644 index 561742815..000000000 --- a/examples/multi-objective/mdp/power/power-timed.nm +++ /dev/null @@ -1,160 +0,0 @@ -// power manager example -mdp - -const int QMAX; // max queue size - -// to model the pm making a choice and then a move being made we need -// two clock ticks for each transition -// first the pm decides tick1 and then the system moves tick2 - -module timer - - c : [0..1]; - - [tick1] c=0 -> (c'=1); - [tick2] c=1 -> (c'=0); - -endmodule - -//------------------------------------------------------------------------- - -// POWER MANAGER -module PM - - pm : [0..4] init 4; - // 0 - go to active - // 1 - go to idle - // 2 - go to idlelp - // 3 - go to stby - // 4 - go to sleep - - [tick1] true -> (pm'=0); - [tick1] true -> (pm'=1); - [tick1] true -> (pm'=2); - [tick1] true -> (pm'=3); - [tick1] true -> (pm'=4); - -endmodule - - -//------------------------------------------------------------------------- - -// SERVICE REQUESTER -module SR - - sr : [0..1] init 0; - // 0 idle - // 1 1req - - [tick2] sr=0 -> 0.898: (sr'=0) + 0.102: (sr'=1); - [tick2] sr=1 -> 0.454: (sr'=0) + 0.546: (sr'=1); - -endmodule - -//------------------------------------------------------------------------- - -// SERVICE PROVIDER - -module SP - - sp : [0..10] init 9; - // 0 active - // 1 idle - // 2 active_idlelp - // 3 idlelp - // 4 idlelp_active - // 5 active_stby - // 6 stby - // 7 stby_active - // 8 active_sleep - // 9 sleep - // 10 sleep_active - - // states where PM has no control (transient states) - [tick2] sp=2 -> 0.75 : (sp'=2) + 0.25 : (sp'=3); // active_idlelp - [tick2] sp=4 -> 0.25 : (sp'=0) + 0.75 : (sp'=4); // idlelp_active - [tick2] sp=5 -> 0.995 : (sp'=5) + 0.005 : (sp'=6); // active_stby - [tick2] sp=7 -> 0.005 : (sp'=0) + 0.995 : (sp'=7); // stby_active - [tick2] sp=8 -> 0.9983 : (sp'=8) + 0.0017 : (sp'=9); // active_sleep - [tick2] sp=10 -> 0.0017 : (sp'=0) + 0.9983 : (sp'=10); // sleep_active - - // states where PM has control - // goto_active - [tick2] sp=0 & pm=0 -> (sp'=0); // active - [tick2] sp=1 & pm=0 -> (sp'=0); // idle - [tick2] sp=3 & pm=0 -> (sp'=4); // idlelp - [tick2] sp=6 & pm=0 -> (sp'=7); // stby - [tick2] sp=9 & pm=0 -> (sp'=10); // sleep - // goto_idle - [tick2] sp=0 & pm=1 -> (sp'=1); // active - [tick2] sp=1 & pm=1 -> (sp'=1); // idle - [tick2] sp=3 & pm=1 -> (sp'=3); // idlelp - [tick2] sp=6 & pm=1 -> (sp'=6); // stby - [tick2] sp=9 & pm=1 -> (sp'=9); // sleep - // goto_idlelp - [tick2] sp=0 & pm=2 -> (sp'=2); // active - [tick2] sp=1 & pm=2 -> (sp'=2); // idle - [tick2] sp=3 & pm=2 -> (sp'=3); // idlelp - [tick2] sp=6 & pm=2 -> (sp'=6); // stby - [tick2] sp=9 & pm=2 -> (sp'=9); // sleep - // goto_stby - [tick2] sp=0 & pm=3 -> (sp'=5); // active - [tick2] sp=1 & pm=3 -> (sp'=5); // idle - [tick2] sp=3 & pm=3 -> (sp'=5); // idlelp - [tick2] sp=6 & pm=3 -> (sp'=6); // stby - [tick2] sp=9 & pm=3 -> (sp'=9); // sleep - // goto_sleep - [tick2] sp=0 & pm=4 -> (sp'=8); // active - [tick2] sp=1 & pm=4 -> (sp'=8); // idle - [tick2] sp=3 & pm=4 -> (sp'=8); // idlelp - [tick2] sp=6 & pm=4 -> (sp'=8); // stby - [tick2] sp=9 & pm=4 -> (sp'=9); // sleep - -endmodule - - -//------------------------------------------------------------------------- - -// SQ -module SQ - - q : [0..QMAX] init 0; - - // serve if busy - [tick2] sr=0 & sp=0 -> (q'=max(q-1,0)); - [tick2] sr=1 & sp=0 -> (q'=q); - - // otherwise do nothing - [tick2] sr=0 & sp>0 -> (q'=q); - [tick2] sr=1 & sp>0 -> (q'=min(q+1,QMAX)); - -endmodule - -//------------------------------------------------------------------------- -//rewards "time" -// [tick2] bat=1 : 1; -//endrewards - -rewards "power" - [tick2] sp=0 & c=1 : 2.5; - [tick2] sp=1 & c=1 : 1.5; - [tick2] sp=2 & c=1 : 2.5; - [tick2] sp=3 & c=1 : 0.8; - [tick2] sp=4 & c=1 : 2.5; - [tick2] sp=5 & c=1 : 2.5; - [tick2] sp=6 & c=1 : 0.3; - [tick2] sp=7 & c=1 : 2.5; - [tick2] sp=8 & c=1 : 2.5; - [tick2] sp=9 & c=1 : 0.1; - [tick2] sp=10 & c=1 : 2.5; -endrewards - -// is an instantaneous property but I suppose we can look at average size -// i.e. divide by the expected number of time steps -rewards "queue" - [tick2] c=1 : q; -endrewards - -rewards "lost" - [tick2] sr=1 & sp>0 & q=2 : 1; -endrewards diff --git a/examples/multi-objective/mdp/power/power-timed.pctl b/examples/multi-objective/mdp/power/power-timed.pctl deleted file mode 100644 index 153141e7c..000000000 --- a/examples/multi-objective/mdp/power/power-timed.pctl +++ /dev/null @@ -1,11 +0,0 @@ -// Average queue size -const double Q; - -// Time bound -const int k; - -// Minimum energy usage over k time-steps, such that average queue size remains below Q -"num_energy": multi(R{"power"}min=? [ C<=k ], R{"queue"}<=Q*k [ C<=k ]) - -// Pareto query: minimum energy usage vs minimum average queue size -"pareto": multi(R{"power"}min=? [ C<=k ], R{"queue"}min=? [ C<=k ]) diff --git a/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi2_time.nm b/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi2_time.nm deleted file mode 100644 index ad0f4a3b7..000000000 --- a/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi2_time.nm +++ /dev/null @@ -1,169 +0,0 @@ -// IPv4: PTA model with digitial clocks -// multi-objective model of the host -// gxn/dxp 28/09/09 - -mdp - -//------------------------------------------------------------- -// VARIABLES -const int N=20; // number of abstract hosts -const int K=2; // number of probes to send - -// PROBABILITIES -const double old = N/65024; // probability pick an ip address being used -//const double old = 1/2; // probability pick an ip address being used -const double new = (1-old); // probability pick a new ip address - -// TIMING CONSTANTS -const int CONSEC = 2; // time interval between sending consecutive probles -const int TRANSTIME = 1; // upper bound on transmission time delay -const int LONGWAIT = 60; // minimum time delay after a high number of address collisions -const int DEFEND = 10; - -const int TIME_MAX_X = 60; // max value of clock x -const int TIME_MAX_Y = 10; // max value of clock y -const int TIME_MAX_Z = 1; // max value of clock z - -// OTHER CONSTANTS -const int MAXCOLL = 10; // maximum number of collisions before long wait - - -//------------------------------------------------------------- -// CONCRETE HOST -module host0 - - x : [0..TIME_MAX_X]; // first clock of the host - y : [0..TIME_MAX_Y]; // second clock of the host - - coll : [0..MAXCOLL]; // number of address collisions - probes : [0..K]; // counter (number of probes sent) - mess : [0..1]; // need to send a message or not - defend : [0..1]; // defend (if =1, try to defend IP address) - - ip : [1..2]; // ip address (1 - in use & 2 - fresh) - - l : [0..4] init 1; // location - // 0 : RECONFIGURE - // 1 : RANDOM - // 2 : WAITSP - // 3 : WAITSG - // 4 : USE - - // RECONFIGURE - [reset] l=0 -> (l'=1); - - // RANDOM (choose IP address) - [rec0] (l=1) -> true; // get message (ignore since have no ip address) - [rec1] (l=1) -> true; // get message (ignore since have no ip address) - // small number of collisions (choose straight away) - [] l=1 & coll 1/3*old : (l'=2) & (ip'=1) & (x'=0) - + 1/3*old : (l'=2) & (ip'=1) & (x'=1) - + 1/3*old : (l'=2) & (ip'=1) & (x'=2) - + 1/3*new : (l'=2) & (ip'=2) & (x'=0) - + 1/3*new : (l'=2) & (ip'=2) & (x'=1) - + 1/3*new : (l'=2) & (ip'=2) & (x'=2); - // large number of collisions: (wait for LONGWAIT) - [time] l=1 & coll=MAXCOLL & x (x'=min(x+1,TIME_MAX_X)); - [] l=1 & coll=MAXCOLL & x=LONGWAIT -> 1/3*old : (l'=2) & (ip'=1) & (x'=0) - + 1/3*old : (l'=2) & (ip'=1) & (x'=1) - + 1/3*old : (l'=2) & (ip'=1) & (x'=2) - + 1/3*new : (l'=2) & (ip'=2) & (x'=0) - + 1/3*new : (l'=2) & (ip'=2) & (x'=1) - + 1/3*new : (l'=2) & (ip'=2) & (x'=2); - - // WAITSP - // let time pass - [time] l=2 & x<2 -> (x'=min(x+1,2)); - // send probe - [send1] l=2 & ip=1 & x=2 & probes (x'=0) & (probes'=probes+1); - [send2] l=2 & ip=2 & x=2 & probes (x'=0) & (probes'=probes+1); - // sent K probes and waited 2 seconds - [configured] l=2 & x=2 & probes=K -> (l'=3) & (probes'=0) & (coll'=0) & (x'=0); - // get message and ip does not match: ignore - [rec0] l=2 & ip!=0 -> (l'=l); - [rec1] l=2 & ip!=1 -> (l'=l); - // get a message with matching ip: reconfigure - [rec1] l=2 & ip=1 -> (l'=0) & (coll'=min(coll+1,MAXCOLL)) & (x'=0) & (probes'=0); - - // WAITSG (sends two gratuitious arp probes) - // time passage - [time] l=3 & mess=0 & defend=0 & x (x'=min(x+1,TIME_MAX_X)); - [time] l=3 & mess=0 & defend=1 & x (x'=min(x+1,TIME_MAX_X)) & (y'=min(y+1,DEFEND)); - - // receive message and same ip: defend - [rec1] l=3 & mess=0 & ip=1 & (defend=0 | y>=DEFEND) -> (defend'=1) & (mess'=1) & (y'=0); - // receive message and same ip: defer - [rec1] l=3 & mess=0 & ip=1 & (defend=0 | y (l'=0) & (probes'=0) & (defend'=0) & (x'=0) & (y'=0); - // receive message and different ip - [rec0] l=3 & mess=0 & ip!=0 -> (l'=l); - [rec1] l=3 & mess=0 & ip!=1 -> (l'=l); - - - // send probe reply or message for defence - [send1] l=3 & ip=1 & mess=1 -> (mess'=0); - [send2] l=3 & ip=2 & mess=1 -> (mess'=0); - // send first gratuitous arp message - [send1] l=3 & ip=1 & mess=0 & x=CONSEC & probes<1 -> (x'=0) & (probes'=probes+1); - [send2] l=3 & ip=2 & mess=0 & x=CONSEC & probes<1 -> (x'=0) & (probes'=probes+1); - // send second gratuitous arp message (move to use) - [send1] l=3 & ip=1 & mess=0 & x=CONSEC & probes=1 -> (l'=4) & (x'=0) & (y'=0) & (probes'=0); - [send2] l=3 & ip=2 & mess=0 & x=CONSEC & probes=1 -> (l'=4) & (x'=0) & (y'=0) & (probes'=0); - - // USE (only interested in reaching this state so do not need to add anything here) - [] l=4 -> true; - -endmodule - -//------------------------------------------------------------- -// error automaton for the environment assumption -// do not get a reply when K probes are sent -const int M; // time between sending and receiving a message - -module env_error2 - - env : [0..1]; // 0 active and 1 done - k : [0..2]; // counts the number of messages sent - c1 : [0..M+1]; // time since first message - c2 : [0..M+1]; // time since second message - error : [0..1]; - - // message with new ip address arrives so done - [send2] error=0 & env=0 -> (env'=1); - // message with old ip address arrives so count - [send1] error=0 & env=0 -> (k'=min(k+1,K)); - // time passgae so update relevant clocks - [time] error=0 & env=0 & k=0 -> true; - [time] error=0 & env=0 & k=1 & min(c1,c2) (c1'=min(c1+1,M+1)); - [time] error=0 & env=0 & k=2 & min(c1,c2) (c1'=min(c1+1,M+1)) & (c2'=min(c2+1,M+1)); - // all clocks reached their bound so an error - [time] error=0 & env=0 & min(c1,c2)=M -> (error'=1); - // send a reply (then done) - [rec1] error=0 & env=0 & k>0 & min(c1,c2)<=M -> (env'=1); - // finished so any action can be performed - [time] error=1 | env=1 -> true; - [send1] error=1 | env=1 -> true; - [send2] error=1 | env=1 -> true; - [rec1] error=1 | env=1 -> true; - -endmodule - -//------------------------------------------------------------- -// error automaton for the time bounded assumption -// host does not send configured signal within T seconds -const int T; - -module time_error - - time_error : [0..1]; - done : [0..1]; - t : [0..T]; - - [time] t (t'=t+1); // time passes and bound not reached - [time] t=T-1 & done=0 & time_error=0 -> (time_error'=1); // bound reached so error - [configured] time_error=0 -> (done'=1); // configured within the time bound - - // when in error or done state can loop with either action - [configured] time_error=1 | done=1 -> true; - [time] time_error=1 | done=1 -> true; - -endmodule diff --git a/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi4_time.nm b/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi4_time.nm deleted file mode 100644 index 71984c286..000000000 --- a/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi4_time.nm +++ /dev/null @@ -1,174 +0,0 @@ -// IPv4: PTA model with digitial clocks -// multi-objective model of the host -// gxn/dxp 28/09/09 - -mdp - -//------------------------------------------------------------- -// VARIABLES -const int N=20; // number of abstract hosts -const int K=4; // number of probes to send - -// PROBABILITIES -const double old = N/65024; // probability pick an ip address being used -//const double old = 0.5; // probability pick an ip address being used -const double new = (1-old); // probability pick a new ip address - -// TIMING CONSTANTS -const int CONSEC = 2; // time interval between sending consecutive probles -const int TRANSTIME = 1; // upper bound on transmission time delay -const int LONGWAIT = 60; // minimum time delay after a high number of address collisions -const int DEFEND = 10; - -const int TIME_MAX_X = 60; // max value of clock x -const int TIME_MAX_Y = 10; // max value of clock y -const int TIME_MAX_Z = 1; // max value of clock z - -// OTHER CONSTANTS -const int MAXCOLL = 10; // maximum number of collisions before long wait - - -//------------------------------------------------------------- -// CONCRETE HOST -module host0 - - x : [0..TIME_MAX_X]; // first clock of the host - y : [0..TIME_MAX_Y]; // second clock of the host - - coll : [0..MAXCOLL]; // number of address collisions - probes : [0..K]; // counter (number of probes sent) - mess : [0..1]; // need to send a message or not - defend : [0..1]; // defend (if =1, try to defend IP address) - - ip : [1..2]; // ip address (1 - in use & 2 - fresh) - - l : [0..4] init 1; // location - // 0 : RECONFIGURE - // 1 : RANDOM - // 2 : WAITSP - // 3 : WAITSG - // 4 : USE - - // RECONFIGURE - [reset] l=0 -> (l'=1); - - // RANDOM (choose IP address) - [rec0] (l=1) -> true; // get message (ignore since have no ip address) - [rec1] (l=1) -> true; // get message (ignore since have no ip address) - // small number of collisions (choose straight away) - [] l=1 & coll 1/3*old : (l'=2) & (ip'=1) & (x'=0) - + 1/3*old : (l'=2) & (ip'=1) & (x'=1) - + 1/3*old : (l'=2) & (ip'=1) & (x'=2) - + 1/3*new : (l'=2) & (ip'=2) & (x'=0) - + 1/3*new : (l'=2) & (ip'=2) & (x'=1) - + 1/3*new : (l'=2) & (ip'=2) & (x'=2); - // large number of collisions: (wait for LONGWAIT) - [time] l=1 & coll=MAXCOLL & x (x'=min(x+1,TIME_MAX_X)); - [] l=1 & coll=MAXCOLL & x=LONGWAIT -> 1/3*old : (l'=2) & (ip'=1) & (x'=0) - + 1/3*old : (l'=2) & (ip'=1) & (x'=1) - + 1/3*old : (l'=2) & (ip'=1) & (x'=2) - + 1/3*new : (l'=2) & (ip'=2) & (x'=0) - + 1/3*new : (l'=2) & (ip'=2) & (x'=1) - + 1/3*new : (l'=2) & (ip'=2) & (x'=2); - - // WAITSP - // let time pass - [time] l=2 & x<2 -> (x'=min(x+1,2)); - // send probe - [send1] l=2 & ip=1 & x=2 & probes (x'=0) & (probes'=probes+1); - [send2] l=2 & ip=2 & x=2 & probes (x'=0) & (probes'=probes+1); - // sent K probes and waited 2 seconds - [configured] l=2 & x=2 & probes=K -> (l'=3) & (probes'=0) & (coll'=0) & (x'=0); - // get message and ip does not match: ignore - [rec0] l=2 & ip!=0 -> (l'=l); - [rec1] l=2 & ip!=1 -> (l'=l); - // get a message with matching ip: reconfigure - [rec1] l=2 & ip=1 -> (l'=0) & (coll'=min(coll+1,MAXCOLL)) & (x'=0) & (probes'=0); - - // WAITSG (sends two gratuitious arp probes) - // time passage - [time] l=3 & mess=0 & defend=0 & x (x'=min(x+1,TIME_MAX_X)); - [time] l=3 & mess=0 & defend=1 & x (x'=min(x+1,TIME_MAX_X)) & (y'=min(y+1,DEFEND)); - - // receive message and same ip: defend - [rec1] l=3 & mess=0 & ip=1 & (defend=0 | y>=DEFEND) -> (defend'=1) & (mess'=1) & (y'=0); - // receive message and same ip: defer - [rec1] l=3 & mess=0 & ip=1 & (defend=0 | y (l'=0) & (probes'=0) & (defend'=0) & (x'=0) & (y'=0); - // receive message and different ip - [rec0] l=3 & mess=0 & ip!=0 -> (l'=l); - [rec1] l=3 & mess=0 & ip!=1 -> (l'=l); - - - // send probe reply or message for defence - [send1] l=3 & ip=1 & mess=1 -> (mess'=0); - [send2] l=3 & ip=2 & mess=1 -> (mess'=0); - // send first gratuitous arp message - [send1] l=3 & ip=1 & mess=0 & x=CONSEC & probes<1 -> (x'=0) & (probes'=probes+1); - [send2] l=3 & ip=2 & mess=0 & x=CONSEC & probes<1 -> (x'=0) & (probes'=probes+1); - // send second gratuitous arp message (move to use) - [send1] l=3 & ip=1 & mess=0 & x=CONSEC & probes=1 -> (l'=4) & (x'=0) & (y'=0) & (probes'=0); - [send2] l=3 & ip=2 & mess=0 & x=CONSEC & probes=1 -> (l'=4) & (x'=0) & (y'=0) & (probes'=0); - - // USE (only interested in reaching this state so do not need to add anything here) - [] l=4 -> true; - -endmodule - -//------------------------------------------------------------- -// error automaton for the environment assumption -// do not get a reply when K probes are sent -const int M; // time between sending and receiving a message - -module env_error4 - - env : [0..1]; // 0 active and 1 done - k : [0..4]; // counts the number of messages sent - c1 : [0..M+1]; // time since first message - c2 : [0..M+1]; // time since second message - c3 : [0..M+1]; // time since third message - c4 : [0..M+1]; // time since fourth message - error : [0..1]; - - // message with new ip address arrives so done - [send2] error=0 & env=0 -> (env'=1); - // message with old ip address arrives so count - [send1] error=0 & env=0 -> (k'=min(k+1,K)); - // time passgae so update relevant clocks - [time] error=0 & env=0 & k=0 -> true; - [time] error=0 & env=0 & k=1 & min(c1,c2,c3,c4) (c1'=min(c1+1,M+1)); - [time] error=0 & env=0 & k=2 & min(c1,c2,c3,c4) (c1'=min(c1+1,M+1)) & (c2'=min(c2+1,M+1)); - [time] error=0 & env=0 & k=3 & min(c1,c2,c3,c4) (c1'=min(c1+1,M+1)) & (c2'=min(c2+1,M+1)) & (c3'=min(c3+1,M+1)); - [time] error=0 & env=0 & k=4 & min(c1,c2,c3,c4) (c1'=min(c1+1,M+1)) & (c2'=min(c2+1,M+1)) & (c3'=min(c3+1,M+1)) & (c4'=min(c4+1,M+1)); - // all clocks reached their bound so an error - [time] error=0 & env=0 & min(c1,c2,c3,c4)=M -> (error'=1); - // send a reply (then done) - [rec1] error=0 & env=0 & k>0 & min(c1,c2,c3,c4)<=M -> (env'=1); - // finished so any action can be performed - [time] error=1 | env=1 -> true; - [send1] error=1 | env=1 -> true; - [send2] error=1 | env=1 -> true; - [send2] error=1 | env=1 -> true; - [rec1] error=1 | env=1 -> true; - -endmodule - -//------------------------------------------------------------- -// error automaton for the time bounded assumption -// host does not send configured signal within T seconds -const int T; - -module time_error - - time_error : [0..1]; - done : [0..1]; - t : [0..T]; - - [time] t (t'=t+1); // time passes and bound not reached - [time] t=T-1 & done=0 & time_error=0 -> (time_error'=1); // bound reached so error - [configured] time_error=0 -> (done'=1); // configured within the time bound - - // when in error or done state can loop with either action - [configured] time_error=1 | done=1 -> true; - [time] time_error=1 | done=1 -> true; - -endmodule diff --git a/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi_time.pctl b/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi_time.pctl deleted file mode 100644 index 078b34c65..000000000 --- a/examples/multi-objective/mdp/zeroconf_time/zeroconf_host_multi_time.pctl +++ /dev/null @@ -1,13 +0,0 @@ -// Max probability of component violating assumption property (checked separately) -const double p_fail = -K=2 ? 0.19 : -K=4 ? 0.006859000000000001 : -K=6 ? 2.476099000000001E-4 : -K=8 ? 8.938717390000006E-6 : -0; - -// Assume-guarantee check via multi-objective -"num_ag": multi(Pmax=? [ F time_error=1 ] , P>=1-p_fail [ G (error=0) ]) - -// Pareto query for assume-guarantee check -"pareto": multi(Pmax=? [ F time_error=1 ] , Pmax=? [ G (error=0) ]) diff --git a/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp b/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp index d4f7920a9..70289b5da 100644 --- a/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp +++ b/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp @@ -9,9 +9,8 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h" -#include "src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h" +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.h" +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h" #include "src/exceptions/NotImplementedException.h" @@ -35,20 +34,17 @@ namespace storm { template std::unique_ptr SparseMdpMultiObjectiveModelChecker::checkMultiObjectiveFormula(CheckTask const& checkTask) { - helper::SparseMultiObjectiveModelCheckerInformation info = helper::SparseMdpMultiObjectivePreprocessingHelper::preprocess(checkTask.getFormula(), this->getModel()); + auto preprocessedData = helper::SparseMdpMultiObjectivePreprocessingHelper::preprocess(checkTask.getFormula(), this->getModel()); std::cout << std::endl; - std::cout << "Preprocessed Information:" << std::endl; - info.printInformationToStream(std::cout); + std::cout << "Preprocessing done." << std::endl; + preprocessedData.printToStream(std::cout); #ifdef STORM_HAVE_CARL - helper::SparseMultiObjectiveModelCheckerHelper::check(info); + helper::SparseMultiObjectiveHelper::check(preprocessedData); #else - STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Multi objective model checking currently requires carl."); + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Multi objective model checking requires carl."); #endif - std::cout << "Information after helper call: " << std::endl; - info.printInformationToStream(std::cout); - return std::unique_ptr(new ExplicitQualitativeCheckResult()); } diff --git a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h deleted file mode 100644 index cc1ea8925..000000000 --- a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h +++ /dev/null @@ -1,71 +0,0 @@ -#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMDPMULTIOBJECTIVEPREPROCESSINGHELPER_H_ -#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMDPMULTIOBJECTIVEPREPROCESSINGHELPER_H_ - -#include "src/logic/Formulas.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h" - -namespace storm { - namespace modelchecker { - - - namespace helper { - - /* - * Helper Class to invoke the necessary preprocessing for multi objective model checking - */ - template - class SparseMdpMultiObjectivePreprocessingHelper { - public: - typedef typename SparseMdpModelType::ValueType ValueType; - typedef typename SparseMdpModelType::RewardModelType RewardModelType; - typedef SparseMultiObjectiveModelCheckerInformation Information; - - /*! - * Preprocesses the given model w.r.t. the given formulas. - * @param originalModel The considered model - * @param originalFormula the considered formula. The subformulas should only contain one OperatorFormula at top level, i.e., the formula is simple. - */ - static Information preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseMdpModelType const& originalModel); - - private: - - /*! - * Inserts the information regarding the given formula (i.e. objective) into info.objectives - * - * @param formula OperatorFormula representing the objective - * @param the information collected so far - */ - static void addObjective(std::shared_ptr const& formula, Information& info); - - /*! - * Sets the timeBound for the given objective - */ - static void setStepBoundOfObjective(typename Information::ObjectiveInformation& currentObjective); - - /*! - * Sets whether we should consider negated rewards - */ - static void setWhetherNegatedRewardsAreConsidered(Information& info); - - /*! - * Apply the neccessary preprocessing for the given formula. - * @param formula the current (sub)formula - * @param info the information collected so far - * @param currentObjective the currently considered objective. The given formula should be a a (sub)formula of this objective - * @param optionalRewardModelName the reward model name that is considered for the formula (if available) - */ - static void preprocessFormula(storm::logic::ProbabilityOperatorFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective); - static void preprocessFormula(storm::logic::RewardOperatorFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective); - static void preprocessFormula(storm::logic::UntilFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective); - static void preprocessFormula(storm::logic::BoundedUntilFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective); - static void preprocessFormula(storm::logic::GloballyFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective); - static void preprocessFormula(storm::logic::EventuallyFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName = boost::none); - static void preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName = boost::none); - - }; - - } - } -} - -#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMDPMULTIOBJECTIVEPREPROCESSINGHELPER_H_ */ diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp new file mode 100644 index 000000000..f90cc6d01 --- /dev/null +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp @@ -0,0 +1,327 @@ +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h" + +#include "src/adapters/CarlAdapter.h" +#include "src/models/sparse/Mdp.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/utility/constants.h" +#include "src/utility/vector.h" +#include "src/settings//SettingsManager.h" +#include "src/settings/modules/MultiObjectiveSettings.h" + +#include "src/exceptions/UnexpectedException.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + + template + typename SparseMultiObjectiveHelper::ReturnType SparseMultiObjectiveHelper::check(PreprocessorData const& data) { + ReturnType result; + result.overApproximation = storm::storage::geometry::Polytope::createUniversalPolytope(); + result.underApproximation = storm::storage::geometry::Polytope::createEmptyPolytope(); + + uint_fast64_t numOfObjectivesWithoutThreshold = 0; + for(auto const& obj : data.objectives) { + if(!obj.threshold) { + ++numOfObjectivesWithoutThreshold; + } + } + if(numOfObjectivesWithoutThreshold == 0) { + achievabilityQuery(data, result); + } else if (numOfObjectivesWithoutThreshold == 1) { + numericalQuery(data, result); + } else if (numOfObjectivesWithoutThreshold == data.objectives.size()) { + paretoQuery(data, result); + } else { + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The number of objecties without threshold is not valid. It should be either 0 (achievabilityQuery), 1 (numericalQuery), or " << data.objectives.size() << " (paretoQuery). Got " << numOfObjectivesWithoutThreshold << " instead."); + } + return result; + } + + template + void SparseMultiObjectiveHelper::achievabilityQuery(PreprocessorData const& data, ReturnType& result) { + Point thresholds; + thresholds.reserve(data.objectives.size()); + storm::storage::BitVector strictThresholds(data.objectives.size()); + for(uint_fast64_t objIndex = 0; objIndex < data.objectives.size(); ++objIndex) { + thresholds.push_back(storm::utility::convertNumber(*data.objectives[objIndex].threshold)); + strictThresholds.set(objIndex, data.objectives[objIndex].thresholdIsStrict); + } + + storm::storage::BitVector individualObjectivesToBeChecked(data.objectives.size(), true); + bool converged=false; + SparseMultiObjectiveWeightVectorChecker weightVectorChecker(data); + do { + WeightVector separatingVector = findSeparatingVector(thresholds, result.underApproximation, individualObjectivesToBeChecked); + refineResult(separatingVector, weightVectorChecker, result); + // Check for convergence + if(!checkIfThresholdsAreSatisfied(overApproximation, thresholds, strictThresholds)){ + result.thresholdsAreAchievable = false; + converged=true; + } + if(checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds)){ + result.thresholdsAreAchievable = true; + converged=true; + } + converged |= (storm::settings::getModule().isMaxIterationsSet() && result.iterations.size() >= storm::settings::getModule().getMaxIterations()); + } while(!converged); + } + + template + void SparseMultiObjectiveHelper::numericalQuery(PreprocessorData const& data, ReturnType& result) { + Point thresholds; + thresholds.reserve(data.objectives.size()); + storm::storage::BitVector strictThresholds(data.objectives.size()); + std::vector> thresholdConstraints; + thresholdConstraints.reserve(data.objectives.size()-1); + uint_fast64_t optimizingObjIndex = data.objectives.size(); // init with invalid value... + for(uint_fast64_t objIndex = 0; objIndex < data.objectives.size(); ++objIndex) { + if(data.objectives[objIndex].threshold) { + thresholds.push_back(storm::utility::convertNumber(*data.objectives[objIndex].threshold)); + WeightVector normalVector(data.objectives.size(), storm::utility::zero()); + normalVector[objIndex] = -storm::utility::one(); + thresholdConstraints.emplace_back(std::move(normalVector), -thresholds.back()); + strictThresholds.set(objIndex, data.objectives[objIndex].thresholdIsStrict); + } else { + optimizingObjIndex = objIndex; + thresholds.push_back(storm::utility::zero()); + } + } + auto thresholdsAsPolytope = storm::storage::geometry::Polytope::create(thresholdConstraints); + WeightVector directionOfOptimizingObjective(data.objectives.size(), storm::utility::zero()); + directionOfOptimizingObjective[optimizingObjIndex] = storm::utility::one(); + + // Try to find one valid solution + storm::storage::BitVector individualObjectivesToBeChecked(data.objectives.size(), true); + individualObjectivesToBeChecked.set(optimizingObjIndex, false); + SparseMultiObjectiveWeightVectorChecker weightVectorChecker(data); + bool converged=false; + do { + WeightVector separatingVector = findSeparatingVector(thresholds, result.underApproximation, individualObjectivesToBeChecked); + refineResult(separatingVector, weightVectorChecker, result); + //Pick the threshold for the optimizing objective low enough so valid solutions are not excluded + thresholds[optimizingObjIndex] = std::min(thresholds[optimizingObjIndex], iterations.back().point[optimizingObjIndex]); + if(!checkIfThresholdsAreSatisfied(overApproximation, thresholds, strictThresholds)){ + result.thresholdsAreAchievable = false; + converged=true; + } + if(checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds)){ + result.thresholdsAreAchievable = true; + converged=true; + } + } while(!converged); + if(*result.thresholdsAreAchievable) { + STORM_LOG_DEBUG("Found a solution that satisfies the objective thresholds."); + individualObjectivesToBeChecked.clear(); + // Improve the found solution. + // Note that we do not have to care whether a threshold is strict anymore, because the resulting optimum should be + // the supremum over all strategies. Hence, one could combine a scheduler inducing the optimum value (but possibly violating strict + // thresholds) and (with very low probability) a scheduler that satisfies all (possibly strict) thresholds. + converged = false; + do { + std::pair optimizationRes = underApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective); + STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "The underapproximation is either unbounded or empty."); + thresholds[optimizingObjIndex] = optimizationRes.first[optimizingObjIndex]; + STORM_LOG_DEBUG("Best solution found so far is " << storm::utility::convertNumber(thresholds[optimizingObjIndex]) << "."); + //Compute an upper bound for the optimum and check for convergence + optimizationRes = overApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective); + if(optimizationRes.second) { + result.precisionOfResult = optimizationRes.first[optimizingObjIndex]; - thresholds[optimizingObjIndex]; + STORM_LOG_DEBUG("Solution can be improved by at most " << storm::utility::convertNumber(*result.precisionOfResult)); + if(storm::utility::convertNumber(*result.precisionOfResult) < storm::settings::getModule().getPrecision()) { + result.numericalResult = thresholds[optimizingObjIndex]; + result.optimumIsAchievable = checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds); + converged=true; + } + } + converged |= (storm::settings::getModule().isMaxIterationsSet() && result.iterations.size() >= storm::settings::getModule().getMaxIterations()); + if(!converged) { + WeightVector separatingVector = findSeparatingVector(thresholds, result.underApproximation, individualObjectivesToBeChecked); + refineResult(separatingVector, weightVectorChecker, result); + } + } while(!converged); + } + } + + template + void SparseMultiObjectiveHelper::paretoQuery(PreprocessorData const& data, ReturnType& result) { + //First consider the objectives individually + for(uint_fast64_t objIndex = 0; objIndex < data.objectives.size(); ++objIndex) { + WeightVector direction(data.objectives.size(), storm::utility::zero()); + direction[objIndex] = storm::utility::one(); + refineResult(direction, weightVectorChecker, result); + } + + bool converged=false; + do { + // Get the halfspace of the underApproximation with maximal distance to a vertex of the overApproximation + std::vector> underApproxHalfspaces = result.underApproximation->getHalfspaces(); + std::vector overApproxVertices = result.overApproximation->getVertices(); + uint_fast64_t farestHalfspaceIndex = underApproxHalfspaces.size(); + RationalNumberType farestDistance = storm::utility::zero(); + for(uint_fast64_t halfspaceIndex = 0; halfspaceIndex < underApproxHalfspaces.size(); ++halfspaceIndex) { + for(auto const& vertex : overApproxVertices) { + RationalNumberType distance = -underApproxHalfspaces[halfspaceIndex].euclideanDistance(vertex); + if(distance > farestDistance) { + farestHalfspaceIndex = halfspaceIndex; + farestDistance = distance; + } + } + } + STORM_LOG_DEBUG("Farest distance between under- and overApproximation is " << storm::utility::convertNumber(farestDistance)); + if(storm::utility::convertNumber(farestDistance) < storm::settings::getModule().getPrecision()) { + result.precisionOfResult = farestDistance; + converged = true; + } + converged |= (storm::settings::getModule().isMaxIterationsSet() && result.iterations.size() >= storm::settings::getModule().getMaxIterations()); + if(!converged) { + WeightVector direction = underApproxHalfspaces[farestHalfspaceIndex].normalVector(); + refineResult(direction, weightVectorChecker, result); + } + } while(!converged); + } + + + template + void SparseMultiObjectiveHelper::refineResult(WeightVector const& direction, SparseWeightedObjectivesModelCheckerHelper const& weightedObjectivesChecker, ReturnType& result) { + weightedObjectivesChecker.check(storm::utility::vector::convertNumericVector(direction)); + STORM_LOG_DEBUG("Result with weighted objectives is " << weightedObjectivesChecker.getInitialStateResultOfObjectives()); + result.iterations.emplace_back(direction, storm::utility::vector::convertNumericVector(weightedObjectivesChecker.getInitialStateResultOfObjectives()), weightedObjectivesChecker.getScheduler()); + + updateOverApproximation(result.iterations.back().point, result.iterations.back().weightVector); + updateUnderApproximation(); + } + + template + typename SparseMultiObjectiveHelper::WeightVector SparseMultiObjectiveHelper::findSeparatingVector(Point const& pointToBeSeparated, std::shared_ptr> const& underApproximation, storm::storage::BitVector& individualObjectivesToBeChecked) const { + STORM_LOG_DEBUG("Searching a weight vector to seperate the point given by " << storm::utility::vector::convertNumericVector(pointToBeSeparated) << "."); + + if(underApproximation->isEmpty()) { + // In this case, every weight vector is separating + uint_fast64_t objIndex = individualObjectivesToBeChecked.getNextSetIndex(0) % pointToBeSeparated.size(); + WeightVector result(pointToBeSeparated.size(), storm::utility::zero()); + result[objIndex] = storm::utility::one(); + individualObjectivesToBeChecked.set(objIndex, false); + return result; + } + + // Reaching this point means that the underApproximation contains halfspaces. + // The seperating vector has to be the normal vector of one of these halfspaces. + // We pick one with maximal distance to the given point. However, weight vectors that correspond to checking individual objectives take precedence. + std::vector> halfspaces = underApproximation->getHalfspaces(); + uint_fast64_t farestHalfspaceIndex = 0; + // Note that we are looking for a halfspace that does not contain the point. Thus, the returned distances are negative. + RationalNumberType farestDistance = -halfspaces[farestHalfspaceIndex].euclideanDistance(pointToBeSeparated); + storm::storage::BitVector nonZeroVectorEntries = ~storm::utility::vector::filterZero(halfspaces[farestHalfspaceIndex].normalVector()); + bool foundSeparatingSingleObjectiveVector = nonZeroVectorEntries.getNumberOfSetBits()==1 && individualObjectivesToBeChecked.get(nonZeroVectorEntries.getNextSetIndex(0)) && farestDistance>storm::utility::zero(); + + for(uint_fast64_t halfspaceIndex = 1; halfspaceIndex < halfspaces.size(); ++halfspaceIndex) { + RationalNumberType distance = -halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated); + nonZeroVectorEntries = ~storm::utility::vector::filterZero(halfspaces[farestHalfspaceIndex].normalVector()); + bool isSeparatingSingleObjectiveVector = nonZeroVectorEntries.getNumberOfSetBits() == 1 && individualObjectivesToBeChecked.get(nonZeroVectorEntries.getNextSetIndex(0)) && distance>storm::utility::zero(); + // Check if this halfspace is 'better' than the current one + if((!foundSeparatingSingleObjectiveVector && isSeparatingSingleObjectiveVector ) || (foundSeparatingSingleObjectiveVector==isSeparatingSingleObjectiveVector && distance>farestDistance)) { + foundSeparatingSingleObjectiveVector |= isSeparatingSingleObjectiveVector; + farestHalfspaceIndex = halfspaceIndex; + farestDistance = distance; + } + } + if(foundSeparatingSingleObjectiveVector) { + nonZeroVectorEntries = ~storm::utility::vector::filterZero(halfspaces[farestHalfspaceIndex].normalVector()); + individualObjectivesToBeChecked.set(nonZeroVectorEntries.getNextSetIndex(0), false); + } + + STORM_LOG_THROW(farestDistance>=storm::utility::zero(), storm::exceptions::UnexpectedException, "There is no seperating vector."); + STORM_LOG_DEBUG("Found separating weight vector: " << storm::utility::vector::convertNumericVector(halfspaces[farestHalfspaceIndex].normalVector()) << "."); + return halfspaces[farestHalfspaceIndex].normalVector(); + } + + template + void SparseMultiObjectiveHelper::updateOverApproximation(std::vector const& iterations, std::shared_ptr>& overApproximation) { + storm::storage::geometry::Halfspace h(iterations.back().weightVector, storm::utility::vector::dotProduct(iterations.back().weightVector, iterations.back().point)); + + // Due to numerical issues, it might be the case that the updated overapproximation does not contain the underapproximation, + // e.g., when the new point is strictly contained in the underapproximation. Check if this is the case. + RationalNumberType maximumOffset = h.offset(); + for(auto const& iteration : iterations){ + maximumOffset = std::max(maximumOffset, storm::utility::vector::dotProduct(h.normalVector(), iteration.point)); + } + if(maximumOffset > h.offset()){ + // We correct the issue by shifting the halfspace such that it contains the underapproximation + h.offset() = maximumOffset; + STORM_LOG_WARN("Numerical issues: The overapproximation would not contain the underapproximation. Hence, a halfspace is shifted by " << storm::utility::convertNumber(h.euclideanDistance(iterations.back().point)) << "."); + } + overApproximation = overApproximation->intersection(h); + STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true)); + } + + template + void SparseMultiObjectiveHelper::updateUnderApproximation(std::vector const& iterations, std::shared_ptr>& underApproximation) { + std::vector paretoPointsVec; + paretoPointsVec.reserve(iterations.size()); + for(auto const& iteration : iterations) { + paretoPointsVec.push_back(iterations.point); + } + + STORM_LOG_WARN("REMOVE ADDING ADDITIONAL VERTICES AS SOON AS HYPRO WORKS FOR DEGENERATED POLYTOPES"); + if(paretoPointsVec.front().size()==2) { + Point p1 = {-10000, -9999}; + Point p2 = {-9999, -10000}; + paretoPointsVec.push_back(p1); + paretoPointsVec.push_back(p2); + } else { + Point p1 = {-10000, -9999, -9999}; + Point p2 = {-9999, -10000, -9999}; + Point p3 = {-9999, -9999, -10000}; + paretoPointsVec.push_back(p1); + paretoPointsVec.push_back(p2); + paretoPointsVec.push_back(p3); + } + + boost::optional upperBounds; + if(!paretoPointsVec.empty()){ + //Get the pointwise maximum of the pareto points + upperBounds = paretoPointsVec.front(); + for(auto paretoPointIt = paretoPointsVec.begin()+1; paretoPointIt != paretoPointsVec.end(); ++paretoPointIt){ + auto upperBoundIt = upperBounds->begin(); + for(auto const& paretoPointCoordinate : *paretoPointIt){ + if(paretoPointCoordinate>*upperBoundIt){ + *upperBoundIt = paretoPointCoordinate; + } + ++upperBoundIt; + } + } + } + + underApproximation = storm::storage::geometry::Polytope::create(paretoPointsVec)->downwardClosure(upperBounds); + STORM_LOG_DEBUG("Updated UnderApproximation to " << underApproximation->toString(true)); + } + + template + bool SparseMultiObjectiveHelper::checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope, Point const& thresholds, storm::storage::BitVector const& strictThresholds) const { + std::vector> halfspaces = polytope->getHalfspaces(); + for(auto const& h : halfspaces) { + RationalNumberType distance = h.distance(thresholds); + if(distance < storm::utility::zero()) { + return false; + } + if(distance == storm::utility::zero()) { + // In this case, the thresholds point is on the boundary of the polytope. + // Check if this is problematic for the strict thresholds + for(auto strictThreshold : strictThresholds) { + if(h.normalVector()[strictThreshold] > storm::utility::zero()) { + return false; + } + } + } + } + return true; + } + +#ifdef STORM_HAVE_CARL + template class SparseMultiObjectiveHelper, storm::RationalNumber>; +#endif + } + } +} diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h new file mode 100644 index 000000000..f6263ace3 --- /dev/null +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h @@ -0,0 +1,67 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEHELPER_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEHELPER_H_ + +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessorReturnType.h" +#include "src/modelchecker/multiObjective/helper/SparseMultiObjectiveWeightVectorChecker.h" +#include "src/storage/geometry/Polytope.h" +#include "src/storage/TotalScheduler.h" + + +namespace storm { + namespace modelchecker { + namespace helper { + + template + class SparseMultiObjectiveHelper { + public: + typedef SpaseMultiObjectivePreprocessorReturnType PreprocessorData; + typedef SpaseMultiObjectiveHelperReturnType ReturnType; + + typedef std::vector Point; + typedef std::vector WeightVector; + + + static ReturnType check(PreprocessorData const& data); + + private: + + static void achievabilityQuery(PreprocessorData const& data, ReturnType& result); + static void numericalQuery(PreprocessorData const& data, ReturnType& result); + static void paretoQuery(PreprocessorData const& data, ReturnType& result); + + /* + * Refines the current result w.r.t. the given direction vector + */ + static void refineResult(WeightVector const& direction, SparseWeightedObjectivesModelCheckerHelper const& weightedObjectivesChecker, ReturnType& result); + + /* + * Returns a weight vector w that separates the under approximation from the given point p, i.e., + * For each x in the under approximation, it holds that w*x <= w*p + * + * @param pointToBeSeparated the point that is to be seperated + * @param underapproximation the current under approximation + * @param objectivesToBeCheckedIndividually stores for each objective whether it still makes sense to check for this objective individually (i.e., with weight vector given by w_{objIndex} = 1 ) + */ + static WeightVector findSeparatingVector(Point const& pointToBeSeparated, std::shared_ptr> const& underApproximation, storm::storage::BitVector& individualObjectivesToBeChecked) const; + + /* + * Updates the over- and under approximation, respectively + * + * @param iterations the iterations performed so far. iterations.back() should be the newest iteration, whose information is not yet included in the approximation. + * @param over/underApproximation the current over/underApproximation that will be updated + */ + static void updateOverApproximation(std::vector const& iterations, std::shared_ptr>& overApproximation); + static void updateUnderApproximation(std::vector const& iterations, std::shared_ptr>& underApproximation); + + /* + * Returns true iff there is a point in the given polytope that satisfies the given thresholds. + * It is assumed that the given polytope contains the downward closure of its vertices. + */ + static bool checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope, Point const& thresholds, storm::storage::BitVector const& strictThresholds) const; + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEHELPER_H_ */ diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelperReturnType.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelperReturnType.h new file mode 100644 index 000000000..2c6e18d40 --- /dev/null +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelperReturnType.h @@ -0,0 +1,56 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEHELPERRETURNTYPE_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEHELPERRETURNTYPE_H_ + +#include +#include +#include +#include + +#include "src/storage/geometry/Polytope.h" +#include "src/storage/TotalScheduler.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + template + struct SparseMultiObjectiveHelperReturnType { + + struct Iteration { + Iteration (std::vector const& weightVector, std::vector const& point, storm::storage::TotalScheduler const& scheduler) : weightVector(weightVector), point(point), scheduler(scheduler) { + //Intentionally left empty + } + + Iteration (std::vector&& weightVector, std::vector&& point, storm::storage::TotalScheduler&& scheduler) : weightVector(weightVector), point(point), scheduler(scheduler) { + //Intentionally left empty + } + + std::vector weightVector; + std::vector point; + storm::storage::TotalScheduler scheduler; + } + + //Stores the results for the individual iterations + std::vector iterations; + + //Stores an over/ under approximation of the set of achievable values + std::shared_ptr> overApproximation; + std::shared_ptr> underApproximation; + + // Stores the result of an achievability query (if applicable) + // For a numerical query, stores whether there is one feasible solution + boost::optional thresholdsAreAchievable; + + //Stores the result of a numerical query (if applicable) + boost::optional numericalResult; + boost::optional optimumIsAchievable; //True if there is an actual scheduler that induces the computed supremum (i.e., supremum == maximum) + + //Stores the achieved precision for numerical and pareto queries + boost::optional precisionOfResult; + + }; + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEHELPERRETURNTYPE_H_ */ diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp deleted file mode 100644 index 677c2fb11..000000000 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp +++ /dev/null @@ -1,490 +0,0 @@ -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h" - -#include "src/adapters/CarlAdapter.h" -#include "src/models/sparse/Mdp.h" -#include "src/models/sparse/StandardRewardModel.h" -#include "src/utility/constants.h" -#include "src/utility/vector.h" -#include "src/settings//SettingsManager.h" -#include "src/settings/modules/MultiObjectiveSettings.h" - - -#include "src/exceptions/UnexpectedException.h" - -namespace storm { - namespace modelchecker { - namespace helper { - - template - SparseMultiObjectiveModelCheckerHelper::SparseMultiObjectiveModelCheckerHelper(Information& info) : info(info), weightedObjectivesChecker(info), numIterations(0) { - overApproximation = storm::storage::geometry::Polytope::createUniversalPolytope(); - underApproximation = storm::storage::geometry::Polytope::createEmptyPolytope(); - } - - template - SparseMultiObjectiveModelCheckerHelper::~SparseMultiObjectiveModelCheckerHelper() { - // Intentionally left empty - } - - template - void SparseMultiObjectiveModelCheckerHelper::check(Information & info) { - SparseMultiObjectiveModelCheckerHelper helper(info); - uint_fast64_t numOfObjectivesWithoutThreshold = 0; - for(auto const& obj : info.objectives) { - if(!obj.threshold) { - ++numOfObjectivesWithoutThreshold; - } - } - if(numOfObjectivesWithoutThreshold == 0) { - helper.achievabilityQuery(); - } else if (numOfObjectivesWithoutThreshold == 1) { - helper.numericalQuery(); - } else if (numOfObjectivesWithoutThreshold == info.objectives.size()) { - helper.paretoQuery(); - } else { - STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The number of objecties without threshold is not valid. It should be either 0 (achievabilityQuery), 1 (numericalQuery), or " << info.objectives.size() << " (paretoQuery). Got " << numOfObjectivesWithoutThreshold << " instead."); - } - } - - template - void SparseMultiObjectiveModelCheckerHelper::achievabilityQuery() { - Point thresholds; - thresholds.reserve(info.objectives.size()); - storm::storage::BitVector strictThresholds(info.objectives.size()); - for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { - STORM_LOG_THROW(info.objectives[objIndex].threshold, storm::exceptions::UnexpectedException, "Can not perform achievabilityQuery: No threshold given for at least one objective."); - thresholds.push_back(storm::utility::convertNumber(*info.objectives[objIndex].threshold)); - strictThresholds.set(objIndex, info.objectives[objIndex].thresholdIsStrict); - } - - storm::storage::BitVector individualObjectivesToBeChecked(info.objectives.size(), true); - bool converged=false; - do { - WeightVector separatingVector = findSeparatingVector(thresholds, individualObjectivesToBeChecked); - refineSolution(separatingVector); - ++numIterations; - // Check for convergence - if(!checkIfThresholdsAreSatisfied(overApproximation, thresholds, strictThresholds)){ - std::cout << "PROPERTY VIOLATED" << std::endl; - converged=true; - } - if(checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds)){ - std::cout << "PROPERTY SATISFIED" << std::endl; - converged=true; - } - converged |= (storm::settings::getModule().isMaxIterationsSet() && numIterations >= storm::settings::getModule().getMaxIterations()); - } while(!converged); - } - - template - void SparseMultiObjectiveModelCheckerHelper::numericalQuery() { - Point thresholds; - thresholds.reserve(info.objectives.size()); - storm::storage::BitVector strictThresholds(info.objectives.size()); - std::vector> thresholdConstraints; - thresholdConstraints.reserve(info.objectives.size()-1); - uint_fast64_t optimizingObjIndex = info.objectives.size(); - for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { - if(info.objectives[objIndex].threshold) { - thresholds.push_back(storm::utility::convertNumber(*info.objectives[objIndex].threshold)); - WeightVector normalVector(info.objectives.size(), storm::utility::zero()); - normalVector[objIndex] = -storm::utility::one(); - thresholdConstraints.emplace_back(std::move(normalVector), -thresholds.back()); - std::cout << "adding threshold constraint" << thresholdConstraints.back().toString(true) << std::endl; - strictThresholds.set(objIndex, info.objectives[objIndex].thresholdIsStrict); - } else { - STORM_LOG_ASSERT(optimizingObjIndex == info.objectives.size(), "Numerical Query invoked but there are multiple objectives without threshold"); - optimizingObjIndex = objIndex; - thresholds.push_back(storm::utility::zero()); - } - } - STORM_LOG_ASSERT(optimizingObjIndex < info.objectives.size(), "Numerical Query invoked but there are no objectives without threshold"); - auto thresholdsAsPolytope = storm::storage::geometry::Polytope::create(thresholdConstraints); - WeightVector directionOfOptimizingObjective(info.objectives.size(), storm::utility::zero()); - directionOfOptimizingObjective[optimizingObjIndex] = storm::utility::one(); - - // Try to find one valid solution - storm::storage::BitVector individualObjectivesToBeChecked(info.objectives.size(), true); - individualObjectivesToBeChecked.set(optimizingObjIndex, false); - do { //TODO introduce convergence criterion? (like max num of iterations) - WeightVector separatingVector = findSeparatingVector(thresholds, individualObjectivesToBeChecked); - refineSolution(separatingVector); - //Pick the threshold for the optimizing objective low enough so valid solutions are not excluded - for(auto const& paretoPoint : paretoOptimalPoints) { - thresholds[optimizingObjIndex] = std::min(thresholds[optimizingObjIndex], paretoPoint.first[optimizingObjIndex]); - } - if(!checkIfThresholdsAreSatisfied(overApproximation, thresholds, strictThresholds)){ - std::cout << "INFEASIBLE OBJECTIVE THRESHOLDS" << std::endl; - return; - } - } while(!checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds)); - STORM_LOG_DEBUG("Found a solution that satisfies the objective thresholds."); - individualObjectivesToBeChecked.clear(); - - // Improve the found solution. - // Note that we do not have to care whether a threshold is strict anymore, because the resulting optimum should be - // the supremum over all strategies. Hence, one could combine a scheduler inducing the optimum value (but possibly violating strict - // thresholds) and (with very low probability) a scheduler that satisfies all (possibly strict) thresholds. - while (true) { - std::pair optimizationRes = underApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective); - STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "The underapproximation is either unbounded or empty."); - thresholds[optimizingObjIndex] = optimizationRes.first[optimizingObjIndex]; - STORM_LOG_DEBUG("Best solution found so far is " << storm::utility::convertNumber(thresholds[optimizingObjIndex]) << "."); - //Compute an upper bound for the optimum and check for convergence - optimizationRes = overApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective); - if(optimizationRes.second) { - RationalNumberType upperBoundOnOptimum = optimizationRes.first[optimizingObjIndex]; - STORM_LOG_DEBUG("Solution can be improved by at most " << storm::utility::convertNumber(upperBoundOnOptimum - thresholds[optimizingObjIndex])); - if(storm::utility::convertNumber(upperBoundOnOptimum - thresholds[optimizingObjIndex]) < storm::settings::getModule().getPrecision()) { std::cout << "Found Optimum: " << storm::utility::convertNumber(thresholds[optimizingObjIndex]) << "." << std::endl; - if(!checkIfThresholdsAreSatisfied(underApproximation, thresholds, strictThresholds)) { - std::cout << "Optimum is only the supremum" << std::endl; - } - return; - } - } - WeightVector separatingVector = findSeparatingVector(thresholds, individualObjectivesToBeChecked); - refineSolution(separatingVector); - } - } - - template - void SparseMultiObjectiveModelCheckerHelper::paretoQuery() { - //First consider the objectives individually - for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { - WeightVector direction(info.objectives.size(), storm::utility::zero()); - direction[objIndex] = storm::utility::one(); - refineSolution(direction); - } - - while(true) { //todo maxNumOfIterations ? - - // Get the halfspace of the underApproximation with maximal distance to a vertex of the overApproximation - std::vector> underApproxHalfspaces = underApproximation->getHalfspaces(); - std::vector overApproxVertices = overApproximation->getVertices(); - uint_fast64_t farestHalfspaceIndex = underApproxHalfspaces.size(); - RationalNumberType farestDistance = storm::utility::zero(); - for(uint_fast64_t halfspaceIndex = 0; halfspaceIndex < underApproxHalfspaces.size(); ++halfspaceIndex) { - for(auto const& vertex : overApproxVertices) { - RationalNumberType distance = underApproxHalfspaces[halfspaceIndex].euclideanDistance(vertex); - // Note that the distances are negative, i.e., we have to consider the minimum distance - if(distance < farestDistance) { - farestHalfspaceIndex = halfspaceIndex; - farestDistance = distance; - } - } - } - STORM_LOG_DEBUG("Farest distance between under- and overApproximation is " << storm::utility::convertNumber(farestDistance)); - if(-storm::utility::convertNumber(farestDistance) < storm::settings::getModule().getPrecision()) { - std::cout << "DONE computing the pareto curve" << std::endl; - return; - } - WeightVector direction = underApproxHalfspaces[farestHalfspaceIndex].normalVector(); - // normalize the direction vector so the entries sum up to one - storm::utility::vector::scaleVectorInPlace(direction, storm::utility::one() / std::accumulate(direction.begin(), direction.end(), storm::utility::zero())); - refineSolution(direction); - } - } - - - template - void SparseMultiObjectiveModelCheckerHelper::refineSolution(WeightVector const& direction) { - weightedObjectivesChecker.check(storm::utility::vector::convertNumericVector(direction)); - - storm::storage::TotalScheduler scheduler = weightedObjectivesChecker.getScheduler(); - Point weightedObjectivesResult = storm::utility::vector::convertNumericVector(weightedObjectivesChecker.getInitialStateResultOfObjectives()); - STORM_LOG_DEBUG("Result with weighted objectives is " << storm::utility::vector::toString(weightedObjectivesResult)); - - // Insert the computed scheduler and check whether we have already seen it before - auto schedulerInsertRes = schedulers.insert(std::make_pair(std::move(scheduler), paretoOptimalPoints.end())); - if(schedulerInsertRes.second){ - // The scheduler is new, so insert the newly computed pareto optimal point. - // For the corresponding scheduler, we safe the entry in the map. - schedulerInsertRes.first->second = paretoOptimalPoints.insert(std::make_pair(std::move(weightedObjectivesResult), std::vector())).first; - } - // In the case where the scheduler is not new, we assume that the corresponding pareto optimal points for the old and new scheduler are equal. - // hence, no additional point will be inserted. - // Note that the values might not be exactly equal due to numerical issues. - - Point const& paretoPoint = schedulerInsertRes.first->second->first; - std::vector& weightVectors = schedulerInsertRes.first->second->second; - weightVectors.push_back(direction); - - updateOverApproximation(paretoPoint, weightVectors.back()); - updateUnderApproximation(); - } - - template - typename SparseMultiObjectiveModelCheckerHelper::WeightVector SparseMultiObjectiveModelCheckerHelper::findSeparatingVector(Point const& pointToBeSeparated, storm::storage::BitVector& individualObjectivesToBeChecked) const { - STORM_LOG_DEBUG("Searching a weight vector to seperate the point given by " << storm::utility::vector::convertNumericVector(pointToBeSeparated) << "."); - // First, we check 'simple' weight vectors that correspond to checking a single objective - while (!individualObjectivesToBeChecked.empty()){ - uint_fast64_t objectiveIndex = individualObjectivesToBeChecked.getNextSetIndex(0); - individualObjectivesToBeChecked.set(objectiveIndex, false); - - WeightVector normalVector(info.objectives.size(), storm::utility::zero()); - normalVector[objectiveIndex] = storm::utility::one(); - - storm::storage::geometry::Halfspace h(normalVector, storm::utility::vector::dotProduct(normalVector, pointToBeSeparated)); - bool hIsSeparating = true; - for(auto const& paretoPoint : paretoOptimalPoints){ - hIsSeparating &= h.contains(paretoPoint.first); - } - if(hIsSeparating) { - return h.normalVector(); - } - } - - if(underApproximation->isEmpty()) { - // In this case [1 0..0] is always separating - WeightVector result(info.objectives.size(), storm::utility::zero()); - result.front() = storm::utility::one(); - return result; - } - - // Get the halfspace of the underApproximation with maximal distance to the given point - // Note that we are looking for a halfspace that does not contain the point. Thus the distances are negative and - // we are actually searching for the one with minimal distance. - std::vector> halfspaces = underApproximation->getHalfspaces(); - uint_fast64_t farestHalfspaceIndex = 0; - RationalNumberType farestDistance = halfspaces[farestHalfspaceIndex].euclideanDistance(pointToBeSeparated); - for(uint_fast64_t halfspaceIndex = 1; halfspaceIndex < halfspaces.size(); ++halfspaceIndex) { - RationalNumberType distance = halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated); - if(distance < farestDistance) { - farestHalfspaceIndex = halfspaceIndex; - farestDistance = distance; - } - } - STORM_LOG_THROW(farestDistance<=storm::utility::zero(), storm::exceptions::UnexpectedException, "There is no seperating halfspace."); - // normalize the result so the entries sum up to one - WeightVector result = halfspaces[farestHalfspaceIndex].normalVector(); - storm::utility::vector::scaleVectorInPlace(result, storm::utility::one() / std::accumulate(result.begin(), result.end(), storm::utility::zero())); - return result; - - - /* old version using LP solving - // We now use LP solving to find a seperating halfspace. - // Note that StoRM's LPSover does not support rational numbers. Hence, we optimize a polytope instead. - // The variables of the LP are the normalVector entries w_0 ... w_{n-1} and the minimal distance d between the halfspace and the paretoOptimalPoints - uint_fast64_t numVariables = pointToBeSeparated.size() + 1; - std::vector> constraints; - constraints.reserve((numVariables-1) + 2 + paretoOptimalPoints.size()); - // w_i >= 0 and d>= 0 <==> -w_i <= 0 and -d <= 0 - for(uint_fast64_t i = 0; i()); - constraintCoefficients[i] = -storm::utility::one(); - RationalNumberType constraintOffset = storm::utility::zero(); - constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset)); - } - // sum_i w_i == 1 <==> sum_i w_i <=1 and sum_i -w_i <= -1 - { - WeightVector constraintCoefficients(numVariables, storm::utility::one()); - constraintCoefficients.back() = storm::utility::zero(); - RationalNumberType constraintOffset = storm::utility::one(); - constraints.emplace_back(constraintCoefficients, constraintOffset); - storm::utility::vector::scaleVectorInPlace(constraintCoefficients, -storm::utility::one()); - constraintOffset = -storm::utility::one(); - constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset)); - } - // let q=pointToBeSeparated. For each paretoPoint x: q*w - x*w >= d <==> (x-q) * w + d <= 0 - for(auto const& paretoPoint : paretoOptimalPoints){ - WeightVector constraintCoefficients(numVariables-1); - storm::utility::vector::subtractVectors(paretoPoint.first, pointToBeSeparated, constraintCoefficients); - constraintCoefficients.push_back(storm::utility::one()); - RationalNumberType constraintOffset = storm::utility::zero(); - constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset)); - } - - WeightVector optimizationVector(numVariables); - optimizationVector.back() = storm::utility::one(); // maximize d - std::pair optimizeResult = storm::storage::geometry::Polytope::create(constraints)->optimize(optimizationVector); - STORM_LOG_THROW(optimizeResult.second, storm::exceptions::UnexpectedException, "There seems to be no seperating halfspace."); - optimizeResult.first.pop_back(); //drop the distance - RationalNumberType offset = storm::utility::vector::dotProduct(optimizeResult.first, pointToBeSeparated); - return storm::storage::geometry::Halfspace(std::move(optimizeResult.first), std::move(offset)); - */ - } - - template - void SparseMultiObjectiveModelCheckerHelper::updateOverApproximation(Point const& newPoint, WeightVector const& newWeightVector) { - storm::storage::geometry::Halfspace h(newWeightVector, storm::utility::vector::dotProduct(newWeightVector, newPoint)); - - // Due to numerical issues, it might be the case that the updated overapproximation does not contain the underapproximation, - // e.g., when the new point is strictly contained in the underapproximation. Check if this is the case. - RationalNumberType maximumOffset = h.offset(); - for(auto const& paretoPoint : paretoOptimalPoints){ - maximumOffset = std::max(maximumOffset, storm::utility::vector::dotProduct(h.normalVector(), paretoPoint.first)); - } - if(maximumOffset > h.offset()){ - // We correct the issue by shifting the halfspace such that it contains the underapproximation - h.offset() = maximumOffset; - STORM_LOG_WARN("Numerical issues: The overapproximation would not contain the underapproximation. Hence, a halfspace is shifted by " << storm::utility::convertNumber(h.euclideanDistance(newPoint)) << "."); - } - overApproximation = overApproximation->intersection(h); - STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true)); - } - - template - void SparseMultiObjectiveModelCheckerHelper::updateUnderApproximation() { - std::vector paretoPointsVec; - paretoPointsVec.reserve(paretoOptimalPoints.size()); - for(auto const& paretoPoint : paretoOptimalPoints) { - paretoPointsVec.push_back(paretoPoint.first); - } - - STORM_LOG_WARN("REMOVE ADDING ADDITIONAL VERTICES AS SOON AS HYPRO WORKS FOR DEGENERATED POLYTOPES"); - Point p1 = {-1000, -999}; - Point p2 = {-999, -1000}; - paretoPointsVec.push_back(p1); - paretoPointsVec.push_back(p2); - - boost::optional upperBounds; - if(!paretoPointsVec.empty()){ - //Get the pointwise maximum of the pareto points - upperBounds = paretoPointsVec.front(); - for(auto paretoPointIt = paretoPointsVec.begin()+1; paretoPointIt != paretoPointsVec.end(); ++paretoPointIt){ - auto upperBoundIt = upperBounds->begin(); - for(auto const& paretoPointCoordinate : *paretoPointIt){ - if(paretoPointCoordinate>*upperBoundIt){ - *upperBoundIt = paretoPointCoordinate; - } - ++upperBoundIt; - } - } - } - - underApproximation = storm::storage::geometry::Polytope::create(paretoPointsVec)->downwardClosure(upperBounds); - STORM_LOG_DEBUG("Updated UnderApproximation to " << underApproximation->toString(true)); - } - - template - bool SparseMultiObjectiveModelCheckerHelper::checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope, Point const& thresholds, storm::storage::BitVector const& strictThresholds) const { - std::vector> halfspaces = polytope->getHalfspaces(); - for(auto const& h : halfspaces) { - RationalNumberType distance = h.distance(thresholds); - if(distance < storm::utility::zero()) { - return false; - } - if(distance == storm::utility::zero()) { - // In this case, the thresholds point is on the boundary of the polytope. - // Check if this is problematic for the strict thresholds - for(auto strictThreshold : strictThresholds) { - if(h.normalVector()[strictThreshold] > storm::utility::zero()) { - return false; - } - } - } - } - return true; - } - - - - /* - - bool existsObjectiveWithoutThreshold = false; - bool existsObjectiveWithStrictThreshold = false; - Point thresholdPoint; - thresholdPoint.reserve(info.objectives.size()); - for(auto const& obj : info.objectives) { - if(obj.threshold) { - thresholdPoint.push_back(storm::utility::convertNumber(*obj.threshold)); - if(obj.thresholdIsStrict) { - existsObjectiveWithStrictThreshold = true; - } - } else { - existsObjectiveWithoutThreshold = true; - thresholdPoint.push_back(storm::utility::one()); - } - } - - if(existsObjectiveWithoutThreshold) { - // We need to find values for the objectives without thresholds in a way that the result is not affected by them. - // We build a polytope that contains valid solutions according to the thresholds and then analyze the intersection with the given polytope. - std::vector> thresholdConstraints; - WeightVector directionForOptimization(info.objectives.size(), storm::utility::zero()); - for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { - if(info.objectives[objIndex].threshold) { - WeightVector normalVector(info.objectives.size(), storm::utility::zero()); - normalVector[objIndex] = -storm::utility::one(); - thresholdConstraints.emplace_back(std::move(normalVector), -thresholdPoint[objIndex]); - directionForOptimization[objIndex] = -storm::utility::one(); - } - } - auto validSolutions = polytope->intersection(storm::storage::geometry::Polytope::create(std::move(thresholdConstraints))); - if(validSolutions->isEmpty()) { - return false; - } - auto optimizationRes = validSolutions->optimize(directionForOptimization); - STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "Couldn't check whether thresholds are satisfied for a given polytope: Didn't expect an infeasible or unbounded solution."); - - for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { - if(!info.objectives[objIndex].threshold) { - thresholdPoint[objIndex] = optimizationRes.first[objIndex] - storm::utility::one(); - } - } - } - - - - - - - // We then check whether the intersection of the thresholdPolytope and the given polytope contains any point (that also satisfies the strict thresholds); - - std::cout << "Valid solutions are " << validSolutions->toString(true) << std::endl; - std::cout << "not empty" << std::endl; - if(existsObjectiveWithStrictThreshold) { - std::cout << "hasObjWithStrictThreshold" << std::endl; - std::cout << "optRes is " << storm::utility::vector::convertNumericVector(optimizationRes.first) << std::endl; - for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { - if(info.objectives[objIndex].threshold && info.objectives[objIndex].thresholdIsStrict) { - RationalNumberType threshold = storm::utility::convertNumber( *info.objectives[objIndex].threshold); - if (optimizationRes.first[objIndex] <= threshold) { - return false; - } - } - } - } - return true; - - - - Point thresholds; - thresholds.reserve(info.objectives.size()); - for(auto const& obj : info.objectives) { - if(obj.threshold) { - thresholds.push_back(storm::utility::convertNumber(*obj.threshold)); - } else { - STORM_LOG_ASSERT(numericalQueryResult, "Checking whether thresholds are satisfied for a numerical query but no threshold for the numerical objective was given."); - thresholds.push_back(*numericalQueryResult); - } - } - std::vector> halfspaces = polytope->getHalfspaces(); - for(auto const& h : halfspaces) { - RationalNumberType distance = h.distance(thresholds); - if(distance < storm::utility::zero()) { - return false; - } - if(distance == storm::utility::zero()) { - // In this case, the thresholds point is on the boundary of the polytope. - // Check if this is problematic for the strict thresholds - for(uint_fast64_t objIndex = 0; objIndex < info.objectives.size(); ++objIndex) { - if(info.objectives[objIndex].thresholdIsStrict && h.normalVector()[objIndex] > storm::utility::zero()) { - return false; - } - } - } - } - return true; - - - - } - */ - -#ifdef STORM_HAVE_CARL - template class SparseMultiObjectiveModelCheckerHelper, storm::RationalNumber>; -#endif - } - } -} diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h deleted file mode 100644 index b8c32ea9d..000000000 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h +++ /dev/null @@ -1,80 +0,0 @@ -#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERHELPER_H_ -#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERHELPER_H_ - -#include -#include -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h" -#include "src/modelchecker/multiObjective/helper/SparseWeightedObjectivesModelCheckerHelper.h" -#include "src/storage/geometry/Polytope.h" -#include "src/storage/TotalScheduler.h" - - -namespace storm { - namespace modelchecker { - namespace helper { - - template - class SparseMultiObjectiveModelCheckerHelper { - public: - typedef typename SparseModelType::ValueType ModelValueType; - typedef typename SparseModelType::RewardModelType RewardModelType; - typedef SparseMultiObjectiveModelCheckerInformation Information; - - - typedef std::vector Point; - typedef std::vector WeightVector; - - static void check(Information& info); - - private: - SparseMultiObjectiveModelCheckerHelper(Information& info); - - ~SparseMultiObjectiveModelCheckerHelper(); - - void achievabilityQuery(); - void numericalQuery(); - void paretoQuery(); - - - - /* - * Refines the solution w.r.t. the given direction vector - */ - void refineSolution(WeightVector const& direction); - - /* - * Returns a weight vector w that separates the underapproximation from the given point p, i.e., - * For each x in the underApproximation, it holds that w*x <= w*p - * - * @param pointToBeSeparated the point that is to be seperated - * @param objectivesToBeCheckedIndividually stores for each objective whether it still makes sense to check for this objective individually (i.e., with weight vector given by w_{objIndex} = 1 ) - */ - WeightVector findSeparatingVector(Point const& pointToBeSeparated, storm::storage::BitVector& individualObjectivesToBeChecked) const; - - void updateOverApproximation(Point const& newPoint, WeightVector const& newWeightVector); - void updateUnderApproximation(); - - /* - * Returns true iff there is a point in the given polytope that satisfies the given thresholds. - * It is assumed that the given polytope contains the downward closure of its vertices. - */ - bool checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope, Point const& thresholds, storm::storage::BitVector const& strictThresholds) const; - - Information& info; - SparseWeightedObjectivesModelCheckerHelper weightedObjectivesChecker; - - //TODO: sort points as needed - std::map> paretoOptimalPoints; - std::unordered_map>::iterator> schedulers; - - uint_fast64_t numIterations; - - std::shared_ptr> overApproximation; - std::shared_ptr> underApproximation; - }; - - } - } -} - -#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERHELPER_H_ */ diff --git a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.cpp similarity index 67% rename from src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp rename to src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.cpp index 36b5a3b42..834977a4d 100644 --- a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.cpp @@ -1,15 +1,11 @@ - #include "src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h" + #include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.h" #include "src/models/sparse/Mdp.h" #include "src/models/sparse/StandardRewardModel.h" - #include "src/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" - -#include "src/transformer/StateDuplicator.h" - #include "src/storage/BitVector.h" - +#include "src/transformer/StateDuplicator.h" #include "src/utility/macros.h" #include "src/utility/vector.h" #include "src/utility/graph.h" @@ -21,24 +17,22 @@ namespace storm { namespace modelchecker { namespace helper { - template - SparseMultiObjectiveModelCheckerInformation SparseMdpMultiObjectivePreprocessingHelper::preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseMdpModelType const& originalModel) { - Information info(std::move(originalModel)); + template + typename SparseMultiObjectivePreprocessor::ReturnType SparseMultiObjectivePreprocessor::preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseModelType const& originalModel) { + ReturnType data(std::move(originalModel)); //Initialize the state mapping. - info.newToOldStateIndexMapping = storm::utility::vector::buildVectorForRange(0, info.preprocessedModel.getNumberOfStates()); + data.newToOldStateIndexMapping = storm::utility::vector::buildVectorForRange(0, data.preprocessedModel.getNumberOfStates()); //Gather information regarding the individual objectives for(auto const& subFormula : originalFormula.getSubFormulas()){ - addObjective(subFormula, info); + addObjective(subFormula, data); } - // Find out whether negated rewards should be considered. - //setWhetherNegatedRewardsAreConsidered(info); //Invoke preprocessing on objectives - for (auto& obj : info.objectives){ + for (auto& obj : data.objectives){ STORM_LOG_DEBUG("Preprocessing objective " << *obj.originalFormula << "."); if(obj.originalFormula->isProbabilityOperatorFormula()){ - preprocessFormula(obj.originalFormula->asProbabilityOperatorFormula(), info, obj); + preprocessFormula(obj.originalFormula->asProbabilityOperatorFormula(), data, obj); } else if(obj.originalFormula->isRewardOperatorFormula()){ - preprocessFormula(obj.originalFormula->asRewardOperatorFormula(), info, obj); + preprocessFormula(obj.originalFormula->asRewardOperatorFormula(), data, obj); } else { STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "Could not preprocess the subformula " << *obj.originalFormula << " of " << originalFormula << " because it is not supported"); } @@ -46,27 +40,27 @@ namespace storm { //We can now remove all original reward models to save some memory std::set origRewardModels = originalFormula.getReferencedRewardModels(); for (auto const& rewModel : origRewardModels){ - info.preprocessedModel.removeRewardModel(rewModel); + data.preprocessedModel.removeRewardModel(rewModel); } - return info; + return data; } - template - void SparseMdpMultiObjectivePreprocessingHelper::addObjective(std::shared_ptr const& formula, Information& info) { + template + void SparseMultiObjectivePreprocessor::addObjective(std::shared_ptr const& formula, ReturnType& data) { STORM_LOG_THROW(formula->isOperatorFormula(), storm::exceptions::InvalidPropertyException, "Expected an OperatorFormula as subformula of multi-objective query but got " << *formula); - typename Information::ObjectiveInformation objective; + ObjectiveInformation objective; objective.originalFormula = formula; - objective.rewardModelName = "objective" + std::to_string(info.objectives.size()); + objective.rewardModelName = "objective" + std::to_string(data.objectives.size()); // Make sure the new reward model gets a unique name - while(info.preprocessedModel.hasRewardModel(objective.rewardModelName)){ + while(data.preprocessedModel.hasRewardModel(objective.rewardModelName)){ objective.rewardModelName = "_" + objective.rewardModelName; } - info.objectives.push_back(std::move(objective)); + data.objectives.push_back(std::move(objective)); } - template - void SparseMdpMultiObjectivePreprocessingHelper::setStepBoundOfObjective(typename Information::ObjectiveInformation& objective){ + template + void SparseMultiObjectivePreprocessor::setStepBoundOfObjective(ObjectiveInformation& objective){ if(objective.originalFormula->isProbabilityOperatorFormula()){ storm::logic::Formula const& f = objective.originalFormula->asProbabilityOperatorFormula().getSubformula(); if(f.isBoundedUntilFormula()){ @@ -84,24 +78,8 @@ namespace storm { } } - template - void SparseMdpMultiObjectivePreprocessingHelper::setWhetherNegatedRewardsAreConsidered(Information& info) { - // Negative Rewards are considered whenever there is an unbounded reward objective that requires to minimize. - // If there is no unbounded reward objective, we make a majority vote. - /* uint_fast64_t numOfMinimizingObjectives = 0; - for(auto const& obj : info.objectives){ - if(obj.originalFormula->isRewardOperatorFormula() && !obj.stepBound){ - info.negatedRewardsConsidered = obj.originalFormulaMinimizes; - return; - } - numOfMinimizingObjectives += obj.originalFormulaMinimizes ? 1 : 0; - } - // Reaching this point means that we make a majority vote - info.negatedRewardsConsidered = (numOfMinimizingObjectives*2) > info.objectives.size(); - */} - - template - void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::ProbabilityOperatorFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective) { + template + void SparseMultiObjectivePreprocessor::preprocessFormula(storm::logic::ProbabilityOperatorFormula const& formula, ReturnType& data, ObjectiveInformation& currentObjective) { //Set the information for the objective if(formula.hasBound()) { currentObjective.threshold = formula.getBound().threshold; @@ -118,26 +96,26 @@ namespace storm { // Invoke preprocessing for subformula if(formula.getSubformula().isUntilFormula()){ - preprocessFormula(formula.getSubformula().asUntilFormula(), info, currentObjective); + preprocessFormula(formula.getSubformula().asUntilFormula(), data, currentObjective); } else if(formula.getSubformula().isBoundedUntilFormula()){ - preprocessFormula(formula.getSubformula().asBoundedUntilFormula(), info, currentObjective); + preprocessFormula(formula.getSubformula().asBoundedUntilFormula(), data, currentObjective); } else if(formula.getSubformula().isGloballyFormula()){ - preprocessFormula(formula.getSubformula().asGloballyFormula(), info, currentObjective); + preprocessFormula(formula.getSubformula().asGloballyFormula(), data, currentObjective); } else if(formula.getSubformula().isEventuallyFormula()){ - preprocessFormula(formula.getSubformula().asEventuallyFormula(), info, currentObjective); + preprocessFormula(formula.getSubformula().asEventuallyFormula(), data, currentObjective); } else { STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The subformula of " << formula << " is not supported."); } } - template - void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::RewardOperatorFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective) { + template + void SparseMultiObjectivePreprocessor::preprocessFormula(storm::logic::RewardOperatorFormula const& formula, ReturnType& data, ObjectiveInformation& currentObjective) { //TODO Make sure that our decision for negative rewards is consistent with the formula - // STORM_LOG_THROW(info.negatedRewardsConsidered == currentObjective.originalFormulaMinimizes, storm::exceptions::InvalidPropertyException, "Decided to consider " << (info.negatedRewardsConsidered ? "negated" : "non-negated") << "rewards, but the formula " << formula << (currentObjective.originalFormulaMinimizes ? " minimizes" : "maximizes") << ". Reward objectives should either all minimize or all maximize."); + // STORM_LOG_THROW(data.negatedRewardsConsidered == currentObjective.originalFormulaMinimizes, storm::exceptions::InvalidPropertyException, "Decided to consider " << (data.negatedRewardsConsidered ? "negated" : "non-negated") << "rewards, but the formula " << formula << (currentObjective.originalFormulaMinimizes ? " minimizes" : "maximizes") << ". Reward objectives should either all minimize or all maximize."); // Check if the reward model is uniquely specified - STORM_LOG_THROW((formula.hasRewardModelName() && info.preprocessedModel.hasRewardModel(formula.getRewardModelName())) - || info.preprocessedModel.hasUniqueRewardModel(), storm::exceptions::InvalidPropertyException, "The reward model is not unique and the formula " << formula << " does not specify a reward model."); + STORM_LOG_THROW((formula.hasRewardModelName() && data.preprocessedModel.hasRewardModel(formula.getRewardModelName())) + || data.preprocessedModel.hasUniqueRewardModel(), storm::exceptions::InvalidPropertyException, "The reward model is not unique and the formula " << formula << " does not specify a reward model."); //Set the information for the objective if(formula.hasBound()) { @@ -155,59 +133,59 @@ namespace storm { // Invoke preprocessing for subformula if(formula.getSubformula().isEventuallyFormula()){ - preprocessFormula(formula.getSubformula().asEventuallyFormula(), info, currentObjective, formula.getOptionalRewardModelName()); + preprocessFormula(formula.getSubformula().asEventuallyFormula(), data, currentObjective, formula.getOptionalRewardModelName()); } else if(formula.getSubformula().isCumulativeRewardFormula()) { - preprocessFormula(formula.getSubformula().asCumulativeRewardFormula(), info, currentObjective, formula.getOptionalRewardModelName()); + preprocessFormula(formula.getSubformula().asCumulativeRewardFormula(), data, currentObjective, formula.getOptionalRewardModelName()); } else { STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The subformula of " << formula << " is not supported."); } } - template - void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::UntilFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective) { + template + void SparseMultiObjectivePreprocessor::preprocessFormula(storm::logic::UntilFormula const& formula, ReturnType& data, ObjectiveInformation& currentObjective) { CheckTask phiTask(formula.getLeftSubformula()); CheckTask psiTask(formula.getRightSubformula()); - storm::modelchecker::SparsePropositionalModelChecker mc(info.preprocessedModel); + storm::modelchecker::SparsePropositionalModelChecker mc(data.preprocessedModel); STORM_LOG_THROW(mc.canHandle(phiTask) && mc.canHandle(psiTask), storm::exceptions::InvalidPropertyException, "The subformulas of " << formula << " should be propositional."); storm::storage::BitVector phiStates = mc.check(phiTask)->asExplicitQualitativeCheckResult().getTruthValuesVector(); storm::storage::BitVector psiStates = mc.check(psiTask)->asExplicitQualitativeCheckResult().getTruthValuesVector(); - auto duplicatorResult = storm::transformer::StateDuplicator::transform(info.preprocessedModel, ~phiStates | psiStates); - info.preprocessedModel = std::move(*duplicatorResult.model); + auto duplicatorResult = storm::transformer::StateDuplicator::transform(data.preprocessedModel, ~phiStates | psiStates); + data.preprocessedModel = std::move(*duplicatorResult.model); // duplicatorResult.newToOldStateIndexMapping now reffers to the indices of the model we had before preprocessing this formula. // This might not be the actual original model. for(auto & originalModelStateIndex : duplicatorResult.newToOldStateIndexMapping){ - originalModelStateIndex = info.newToOldStateIndexMapping[originalModelStateIndex]; + originalModelStateIndex = data.newToOldStateIndexMapping[originalModelStateIndex]; } - info.newToOldStateIndexMapping = std::move(duplicatorResult.newToOldStateIndexMapping); + data.newToOldStateIndexMapping = std::move(duplicatorResult.newToOldStateIndexMapping); // build stateAction reward vector that gives (one*transitionProbability) reward whenever a transition leads from the firstCopy to a psiState - storm::storage::BitVector newPsiStates(info.preprocessedModel.getNumberOfStates(), false); + storm::storage::BitVector newPsiStates(data.preprocessedModel.getNumberOfStates(), false); for(auto const& oldPsiState : psiStates){ //note that psiStates are always located in the second copy newPsiStates.set(duplicatorResult.secondCopyOldToNewStateIndexMapping[oldPsiState], true); } - std::vector objectiveRewards(info.preprocessedModel.getTransitionMatrix().getRowCount(), storm::utility::zero()); + std::vector objectiveRewards(data.preprocessedModel.getTransitionMatrix().getRowCount(), storm::utility::zero()); for (auto const& firstCopyState : duplicatorResult.firstCopy) { - for (uint_fast64_t row = info.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[firstCopyState]; row < info.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[firstCopyState + 1]; ++row) { - objectiveRewards[row] = info.preprocessedModel.getTransitionMatrix().getConstrainedRowSum(row, newPsiStates); + for (uint_fast64_t row = data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[firstCopyState]; row < data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[firstCopyState + 1]; ++row) { + objectiveRewards[row] = data.preprocessedModel.getTransitionMatrix().getConstrainedRowSum(row, newPsiStates); } } if(currentObjective.isNegative){ storm::utility::vector::scaleVectorInPlace(objectiveRewards, -storm::utility::one()); } - info.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); + data.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); } - template - void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::BoundedUntilFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective) { + template + void SparseMultiObjectivePreprocessor::preprocessFormula(storm::logic::BoundedUntilFormula const& formula, ReturnType& data, ObjectiveInformation& currentObjective) { STORM_LOG_THROW(formula.hasDiscreteTimeBound(), storm::exceptions::InvalidPropertyException, "Expected a boundedUntilFormula with a discrete time bound but got " << formula << "."); currentObjective.stepBound = formula.getDiscreteTimeBound(); - preprocessFormula(storm::logic::UntilFormula(formula.getLeftSubformula().asSharedPointer(), formula.getRightSubformula().asSharedPointer()), info, currentObjective); + preprocessFormula(storm::logic::UntilFormula(formula.getLeftSubformula().asSharedPointer(), formula.getRightSubformula().asSharedPointer()), data, currentObjective); } - template - void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::GloballyFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective) { + template + void SparseMultiObjectivePreprocessor::preprocessFormula(storm::logic::GloballyFormula const& formula, ReturnType& data, ObjectiveInformation& currentObjective) { // The formula will be transformed to an until formula for the complementary event currentObjective.isInverted = true; if(currentObjective.threshold) { @@ -221,54 +199,54 @@ namespace storm { } currentObjective.isNegative = !currentObjective.isNegative; auto negatedSubformula = std::make_shared(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, formula.getSubformula().asSharedPointer()); - preprocessFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), negatedSubformula), info, currentObjective); + preprocessFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), negatedSubformula), data, currentObjective); } - template - void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::EventuallyFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName) { + template + void SparseMultiObjectivePreprocessor::preprocessFormula(storm::logic::EventuallyFormula const& formula, ReturnType& data, ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName) { if(formula.isReachabilityProbabilityFormula()){ - preprocessFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), formula.getSubformula().asSharedPointer()), info, currentObjective); + preprocessFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), formula.getSubformula().asSharedPointer()), data, currentObjective); return; } STORM_LOG_THROW(formula.isReachabilityRewardFormula(), storm::exceptions::InvalidPropertyException, "The formula " << formula << " neither considers reachability Probabilities nor reachability rewards"); CheckTask targetTask(formula.getSubformula()); - storm::modelchecker::SparsePropositionalModelChecker mc(info.preprocessedModel); + storm::modelchecker::SparsePropositionalModelChecker mc(data.preprocessedModel); STORM_LOG_THROW(mc.canHandle(targetTask), storm::exceptions::InvalidPropertyException, "The subformula of " << formula << " should be propositional."); storm::storage::BitVector targetStates = mc.check(targetTask)->asExplicitQualitativeCheckResult().getTruthValuesVector(); STORM_LOG_WARN("There is no check for reward finiteness yet"); //TODO - auto duplicatorResult = storm::transformer::StateDuplicator::transform(info.preprocessedModel, targetStates); - info.preprocessedModel = std::move(*duplicatorResult.model); + auto duplicatorResult = storm::transformer::StateDuplicator::transform(data.preprocessedModel, targetStates); + data.preprocessedModel = std::move(*duplicatorResult.model); // duplicatorResult.newToOldStateIndexMapping now reffers to the indices of the model we had before preprocessing this formula. // This might not be the actual original model. for(auto & originalModelStateIndex : duplicatorResult.newToOldStateIndexMapping){ - originalModelStateIndex = info.newToOldStateIndexMapping[originalModelStateIndex]; + originalModelStateIndex = data.newToOldStateIndexMapping[originalModelStateIndex]; } - info.newToOldStateIndexMapping = std::move(duplicatorResult.newToOldStateIndexMapping); + data.newToOldStateIndexMapping = std::move(duplicatorResult.newToOldStateIndexMapping); // Add a reward model that gives zero reward to the states of the second copy. - std::vector objectiveRewards = info.preprocessedModel.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : "").getTotalRewardVector(info.preprocessedModel.getTransitionMatrix()); + std::vector objectiveRewards = data.preprocessedModel.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : "").getTotalRewardVector(data.preprocessedModel.getTransitionMatrix()); storm::utility::vector::setVectorValues(objectiveRewards, duplicatorResult.secondCopy, storm::utility::zero()); if(currentObjective.isNegative){ storm::utility::vector::scaleVectorInPlace(objectiveRewards, -storm::utility::one()); } - info.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); + data.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); } - template - void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName) { + template + void SparseMultiObjectivePreprocessor::preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, ReturnType& data, ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName) { STORM_LOG_THROW(formula.hasDiscreteTimeBound(), storm::exceptions::InvalidPropertyException, "Expected a cumulativeRewardFormula with a discrete time bound but got " << formula << "."); currentObjective.stepBound = formula.getDiscreteTimeBound(); - std::vector objectiveRewards = info.preprocessedModel.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : "").getTotalRewardVector(info.preprocessedModel.getTransitionMatrix()); + std::vector objectiveRewards = data.preprocessedModel.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : "").getTotalRewardVector(data.preprocessedModel.getTransitionMatrix()); if(currentObjective.isNegative){ storm::utility::vector::scaleVectorInPlace(objectiveRewards, -storm::utility::one()); } - info.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); + data.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); } - template class SparseMdpMultiObjectivePreprocessingHelper>; + template class SparseMultiObjectivePreprocessor>; } diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.h new file mode 100644 index 000000000..ee0575c04 --- /dev/null +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.h @@ -0,0 +1,67 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEPREPROCESSOR_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEPREPROCESSOR_H_ + +#include "src/logic/Formulas.h" +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessorReturnType.h" + +namespace storm { + namespace modelchecker { + + + namespace helper { + + /* + * Helper Class to invoke the necessary preprocessing for multi objective model checking + */ + template + class SparseMultiObjectivePreprocessor { + public: + typedef typename SparseModelType::ValueType ValueType; + typedef typename SparseModelType::RewardModelType RewardModelType; + typedef SparseMultiObjectivePreprocessorReturnType ReturnType; + typedef typename SparseMultiObjectivePreprocessorReturnType::ObjectiveInformation ObjectiveInformation; + + /*! + * Preprocesses the given model w.r.t. the given formulas. + * @param originalModel The considered model + * @param originalFormula the considered formula. The subformulas should only contain one OperatorFormula at top level, i.e., the formula is simple. + */ + static ReturnType preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseModelType const& originalModel); + + private: + + /*! + * Inserts the information regarding the given formula (i.e. objective) into info.objectives + * + * @param formula OperatorFormula representing the objective + * @param the information collected so far + */ + static void addObjective(std::shared_ptr const& formula, ReturnType& data); + + /*! + * Sets the timeBound for the given objective + */ + static void setStepBoundOfObjective(ObjectiveInformation& currentObjective); + + /*! + * Apply the neccessary preprocessing for the given formula. + * @param formula the current (sub)formula + * @param data the information collected so far + * @param currentObjective the currently considered objective. The given formula should be a a (sub)formula of this objective + * @param optionalRewardModelName the reward model name that is considered for the formula (if available) + */ + static void preprocessFormula(storm::logic::ProbabilityOperatorFormula const& formula, ReturnType& info, ObjectiveInformation& currentObjective); + static void preprocessFormula(storm::logic::RewardOperatorFormula const& formula, ReturnType& info, ObjectiveInformation& currentObjective); + static void preprocessFormula(storm::logic::UntilFormula const& formula, ReturnType& info, ObjectiveInformation& currentObjective); + static void preprocessFormula(storm::logic::BoundedUntilFormula const& formula, ReturnType& info, ObjectiveInformation& currentObjective); + static void preprocessFormula(storm::logic::GloballyFormula const& formula, ReturnType& info, ObjectiveInformation& currentObjective); + static void preprocessFormula(storm::logic::EventuallyFormula const& formula, ReturnType& info, ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName = boost::none); + static void preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, ReturnType& info, ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName = boost::none); + + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEPREPROCESSOR_H_ */ diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessorReturnType.h similarity index 91% rename from src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h rename to src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessorReturnType.h index fd64a6468..e9e198778 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessorReturnType.h @@ -1,9 +1,10 @@ -#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERINFORMATION_H_ -#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERINFORMATION_H_ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEPREPROCESSORRETURNTYPE_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEPREPROCESSORRETURNTYPE_H_ #include #include #include +#include #include "src/logic/Formulas.h" @@ -28,7 +29,7 @@ namespace storm { bool thresholdIsStrict = false; boost::optional stepBound; - void printInformationToStream(std::ostream& out) const { + void printToStream(std::ostream& out) const { out << std::setw(30) << originalFormula->toString(); out << " \t("; out << (isNegative ? "-x, " : " "); @@ -52,26 +53,19 @@ namespace storm { } }; - SparseMultiObjectiveModelCheckerInformation(SparseModelType const& model) : preprocessedModel(model), originalModel(model) { //Intentionally left empty } - - SparseModelType preprocessedModel; - SparseModelType const& originalModel; - std::vector newToOldStateIndexMapping; - std::vector objectives; - - void printInformationToStream(std::ostream& out) { + void printToStream(std::ostream& out) { out << "---------------------------------------------------------------------------------------------------------------------------------------" << std::endl; - out << " Multi-objective Model Checker Information " << std::endl; + out << " Multi-objective Preprocessor Result " << std::endl; out << "---------------------------------------------------------------------------------------------------------------------------------------" << std::endl; out << std::endl; out << "Objectives:" << std::endl; out << "--------------------------------------------------------------" << std::endl; for (auto const& obj : objectives) { - obj.printInformationToStream(out); + obj.printToStream(out); } out << "--------------------------------------------------------------" << std::endl; out << std::endl; @@ -84,10 +78,14 @@ namespace storm { out << "---------------------------------------------------------------------------------------------------------------------------------------" << std::endl; } + SparseModelType preprocessedModel; + SparseModelType const& originalModel; + std::vector newToOldStateIndexMapping; + std::vector objectives; }; } } } -#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERINFORMATION_H_ */ +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEPREPROCESSORRETURNTYPE_H_ */ diff --git a/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp similarity index 98% rename from src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.cpp rename to src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp index 5f97dbdba..14a20f23c 100644 --- a/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp @@ -144,7 +144,7 @@ namespace storm { //std::cout << "stateActionRewardVector for objective " << objIndex << " is " << storm::utility::vector::toString(info.preprocessedModel.getRewardModel(info.objectives[objIndex].rewardModelName).getStateActionRewardVector()) << std::endl; //std::cout << "deterministic state rewards for objective " << objIndex << " are " << storm::utility::vector::toString(deterministicStateRewards) << std::endl; - storm::storage::BitVector statesWithRewards = storm::utility::vector::filter(deterministicStateRewards, [&] (ValueType const& value) -> bool {return !storm::utility::isZero(value);}); + storm::storage::BitVector statesWithRewards = ~storm::utility::vector::filterZero(deterministicStateRewards); // As target states, we take the states from which no reward is reachable. STORM_LOG_WARN("TODO: target state selection is currently only valid for reachability properties..."); //TODO: we should be able to give some hint to the solver.. diff --git a/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.h similarity index 89% rename from src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.h rename to src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.h index 94f7b9e76..f6b2600bd 100644 --- a/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.h +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.h @@ -1,9 +1,9 @@ -#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEWEIGHTEDOBJECTIVESMODELCHECKERHELPER_H_ -#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEWEIGHTEDOBJECTIVESMODELCHECKERHELPER_H_ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEWEIGHTVECTORCHECKER_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEWEIGHTVECTORCHECKER_H_ #include -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h" +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessorReturnType.h" #include "src/storage/TotalScheduler.h" namespace storm { @@ -11,13 +11,13 @@ namespace storm { namespace helper { /*! - * Helper Class that takes a MultiObjectiveInformation and a weight vector and ... + * Helper Class that takes preprocessed multi objective data and a weight vector and ... * - computes the maximal expected reward w.r.t. the weighted sum of the rewards of the individual objectives * - extracts the scheduler that induces this maximum * - computes for each objective the value induced by this scheduler */ template - class SparseWeightedObjectivesModelCheckerHelper { + class SparseMultiObjectiveWeightVectorChecker { public: typedef typename SparseModelType::ValueType ValueType; typedef typename SparseModelType::RewardModelType RewardModelType; @@ -87,4 +87,4 @@ namespace storm { } } -#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEWEIGHTEDOBJECTIVESMODELCHECKERHELPER_H_ */ +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEWEIGHTEDVECTORCHECKER_H_ */ diff --git a/src/utility/vector.h b/src/utility/vector.h index 37500feb7..4b5a26687 100644 --- a/src/utility/vector.h +++ b/src/utility/vector.h @@ -405,6 +405,18 @@ namespace storm { return filter(values, fnc); } + /*! + * Retrieves a bit vector containing all the indices for which the value at this position is equal to zero + * + * @param values The vector of values. + * @return The resulting bit vector. + */ + template + storm::storage::BitVector filterZero(std::vector const& values) { + std::function fnc = [] (T const& value) -> bool { return storm::utility::isZero(value); }; + return filter(values, fnc); + } + /** * Sum the entries from values that are set to one in the filter vector. * @param values