-- Copyright (C) International Computer Science Institute, 1995.  COPYRIGHT  --
-- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
-- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in    --
-- the file "Doc/License" of the Sather distribution.  The license is also   --
-- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA.  --
--------> Please email comments to "sather-bugs@icsi.berkeley.edu". <----------
-- Implements the DE algorithm in Sather
-- pSather version by Claudio Fleiner
-- Note that the original source code was spread over multiple
-- files, but to simplify retrieval over the WWW those
-- files have been merged
class DE{T<$DE_PROBLEM} is
   private attr D:INT;
   private attr NP:INT;
   private attr population:AREF{T};
   private attr errors:AREF{FLT};
   private attr new_population:AREF{T};
   private attr new_errors:AREF{FLT};
   private attr max_error:FLT;
   private attr id:INT;
   private attr generation:INT;
   private attr rnd:AREF{RND};
   private attr errt:AREF{FLT};
   private attr idt:AREF{INT};
   shared generations:INT:=1000;
   shared print_result:BOOL:=false;
   
   create(vec_size:INT,pop_size:INT):SAME is
      r::=new;
      r.D:=vec_size;
      r.NP:=pop_size;
      r.new_population:=#AREF{T}(pop_size); 
      loop r.new_population[0.upto!(r.NP-1)]:=#T(r.D); end;
      if print_result then #OUT+r.new_population[pop_size-1][vec_size-1]+"\n"; end;
      r.population:=#AREF{T}(pop_size);
      loop r.population[0.upto!(r.NP-1)]:=#T(r.D); end;
      r.errors:=#AREF{FLT}(pop_size);
      r.new_errors:=#AREF{FLT}(pop_size);
      r.id:= -1;
      r.rnd:=#AREF{RND}(THREADS);
      r.errt:=#AREF{FLT}(THREADS);
      r.idt:=#AREF{INT}(THREADS);
      return r;
   end;

   create_random_population(size:FLT,seed:INT) is
      loop
	 rnd.aset!(#RND);
	 -- rnd.aelt!.init(seed-1-rnd.aind!);
      end;

      half::=#POLY_FIT(D);
      half.to(-0.5);
      loop
      	 act::=population.aelt!;
	 j::=population.aind!;
	 loop
	    i::=act.ind!;
	    act[i]:=2.0*size*(RND::uniform.flt-0.5);
	 end;
	 errors.aset!(act.error);
	 if j=0 then
	    max_error:=errors[0];
	    id:=0;
	 else
	    if max_error>errors[j] then
	       max_error:=errors[j];
	       id:=j;
	    end;
	 end;
      end;
      generation:=1;
      if print_result then #OUT+generation+". "+max_error+": "+population[id].str+"\n"; end;
   end;

   print is 
      loop
         act::=population.aelt!;
	 act2::=new_population.aelt!;
--	 #OUT+errors.aelt!+" : "+act.str+"\n"+new_errors.aelt!+" : "+act2.str+"\n";
	 if print_result then #OUT+"--- "+errors.aelt!+" : "+act.str+"\n"; end;
      end;
   end;

   next_generation(t:INT,CR:FLT,F:FLT,from:INT,to:INT) is
      loop
	 r1,r2,r3:INT;
	 i::=from.upto!(to);
	 act::=#T(D);
	 act.to(population[i]);
	 loop r1:=rnd[t].int(0,NP-1);while!(r1=i); end;
	 loop r2:=rnd[t].int(0,NP-1);while!(r2=r1 or r2=i); end;
	 loop r3:=rnd[t].int(0,NP-1);while!(r3=r1 or r3=r2 or r3=i); end;
	 j::=rnd[t].int(0,D-1);
	 loop
	    act[j]:=population[r1][j]+F*(population[r2][j]-population[r3][j]);
	    j:=(j+1)%D;
	    dummy::=act.ind!;
	    while!(rnd[t].uniform.flt<CR);
	 end;
	 err::=act.error;
	 if err<errors[i] then
	    new_population[i].to(act);
	    new_errors[i]:=err;
	    if err<errt[t] then
	       errt[t]:=err;
	       idt[t]:=i;
	    end;
	 else
	    new_population[i].to(population[i]);
	    new_errors[i]:=errors[i];
         end;
      end;
   end;

   consolidate(g:INT) is
      old::=population;population:=new_population;new_population:=old;
      olde::=errors;errors:=new_errors;new_errors:=olde;
      generation:=generation+1;
      loop 
	 i::=errt.aind!;
	 if errt[i]<max_error then
	    max_error:=errt[i];
	    id:=idt[i];
	    if print_result then #OUT+generation+". "+max_error+": "+population[id].str+"\n"; end;
	 else
	    errt[i]:=max_error;
	 end;
      end;
      if g=generations-1 and print_result then #OUT+generation+". "+max_error+": "+population[id].str+"\n"; end;
      -- (#OUT+generation+". "+max_error+"\r").flush;
   end;

   private const THREADS:=cluster_size;

   next_generation(CR:FLT,F:FLT) is
      loop
          errt.aset!(max_error);
      end;
      parloop
         t::=0.upto!(THREADS-1);
      do
	 if print_result then #OUT+("thread "+THREAD_ID::me.str+"\n"); end;
	 loop g::=generations.times!;
	    next_generation(t,CR,F,(NP*t)/THREADS,(NP*(t+1))/THREADS-1);
	    sync;
	    if t=0 then consolidate(g); end;
	    sync;
         end;
      end;
   end;
end;

-- vim:sw=3:nosmartindent
abstract class $DE_PROBLEM is
   error:FLT;
   aset(ind:INT,val:FLT);
   aget(ind:INT):FLT;
end;

class POLY_FIT < $DE_PROBLEM  is
   include VEC ; -- times_add_into->, times->, times_in->, outer_product->;
   shared points:INT:=600;

   eval(a:FLT):FLT is
      r::=0.0;
      loop r:=r*a+elt!; end;
      return r;
   end;

   error:FLT is -- calculate the error function
      r::=0.0;
      c::=0.0;
      incr::=2.0/points.flt;
      x::=-1.0;
      loop i::=0.upto!(points);
         c:=eval(x);
	 if c<-1.0 or c>=1.0 then
	    r:=r+(1.0-c)*(1.0-c);
	 end;
	 x:=x+incr;
      end;
      c:=eval(1.2)-72.661;
      if c<0.0 then 
	 r:=r+c*c;
      end;
      c:=eval(-1.2)-72.661;
      if c<0.0 then 
	 r:=r+c*c;
      end;
      return r;
   end;
end;

class MAIN is
   main(arg:ARRAY{STR}):INT is
      -- GC::concurrency(4);
      pop::=30;
      if arg.size>=2 and arg[1]="-p" then
	 DE{POLY_FIT}::print_result:=true;
	 a::=#ARRAY{STR}(arg.size-1);
	 a[0]:=arg[0];
	 loop a[1.upto!(a.size-1)]:=arg[2.up!]; end;
	 arg:=a;
      end;
      if arg.size=2 then
	 POLY_FIT::points:=#INT(arg[1]);
      elsif arg.size=3 then
	 POLY_FIT::points:=#INT(arg[1]);
	 DE{POLY_FIT}::generations:=#INT(arg[2]);
      elsif arg.size=4 then
	 POLY_FIT::points:=#INT(arg[1]);
	 DE{POLY_FIT}::generations:=#INT(arg[2]);
	 pop:=#INT(arg[3]);
      elsif arg.size/=1 then
         #ERR+"Usage: "+arg[0]+"[accuracy:60-60000,600 [generations:1-10000,1000 [population:10-300,30]]]\n";
	 return 1;
      end;

      p::=#DE{POLY_FIT}(9,pop);
      p.create_random_population(100.0,-3);
      p.print;
      t::=#TIMES;
      p.next_generation(1.0,0.8);
      #OUT+t.elapsed.str;
      return 0;
   end;
end;