## simple MultiStart for DGP
##   assuming realization is in var x{V,K} <= M, >= -M
option randseed 0; #changes the random seed to a number based on the system clock, which will be a different every time the command is given. This changes all random parameters each time
option solver_msg 0;

include "dgp_ms.opt";

### nonconfigurable parameters

# termination 
param ms_termination binary default 0;
# realization error wrt given distances
param ms_err;
# best solution so far
param ms_xstar{V,K};
# lowest realization error so far
param ms_errstar;
# CPU time
param ms_cputime default 0;

# solve the problem locally
let{v in V, k in K} x[v,k] := Uniform(orig_xL[v,k], orig_xU[v,k]);
solve > nul;  # shut off the output competely
# initialize first "best so far" point
let{v in V, k in K} ms_xstar[v,k] := x[v,k];
# compute the realization error
let ms_errstar := 
   (1/card(E))*sum{(u,v) in E} 
      abs(sqrt(abs(sum{k in K} (x[u,k]-x[v,k])^2)) - c[u,v]);
printf "ms: starting solution err = %g\n", ms_errstar;

param ms_opt binary, default 0;
param ms_k integer, default 0;
repeat while(ms_termination == 0) {
  # initialize iteration counter at each loop
  let ms_k := ms_k + 1;
  printf "ms: itn=%d", ms_k;
  # choose initial random point
  let{v in V, k in K} x[v,k] := Uniform(-M, M);
  # solve the problem locally
  solve > nul;
  # compute the realization error
  let ms_err := 
   (1/card(E))*sum{(u,v) in E} 
      abs(sqrt(abs(sum{k in K} (x[u,k]-x[v,k])^2)) - c[u,v]);
  # check for global optimality
  if (ms_err < ms_epsilon) then {
    printf " global optimum found, err = %g", ms_err;
    let {v in V, k in K} ms_xstar[v,k] := x[v,k];
    let ms_errstar := ms_err; 
    let ms_opt := 1;
    let ms_termination := 1;
    break;
  } else if (ms_err < ms_errstar) then {
    printf " improvement by %g at iteration %d,", ms_err, ms_k;
    # error is lower, x is now best so far
    let {v in V, k in K} ms_xstar[v,k] := x[v,k];
    let ms_errstar := ms_err; 
  }   
  # check termination condition
  let ms_cputime := _ampl_user_time + _total_solve_user_time;   # user_CPU_seconds_used_by_the_AMPL_process_itself + user_CPU_seconds_used_by_all_solve_commands
  printf " cputime = %g\n", ms_cputime;
  if (ms_cputime > ms_maxcputime) then {
    let ms_termination := 1; 
  }
  if (ms_termination == 1) then {
    printf "ms: terminating, max cpu time (%g) exceeded\n", ms_maxcputime;
  }
}

# save the solution into x
let {v in V, k in K} x[v,k] := ms_xstar[v,k];