mdecourse / vecp2018

車輛工程系大一計算機程式 (課程已經結束)
GNU Affero General Public License v3.0
0 stars 0 forks source link

全球資訊網上的 Lua 程式應用 #9

Closed mdecourse closed 6 years ago

mdecourse commented 6 years ago

使用 Javascript 改寫 Lua 解譯器之後

請將下列六行程式複製到 http://mde.tw/vecp2018/lua/SchoolProject/ 中執行:

-- 利用 dofile() 從網頁中取得 star() 函式的定義, 並對應到 local 變數 star 中
local star = dofile("https://raw.githubusercontent.com/mdecourse/vecp2018_1/gh-pages/w5/star.lua").star
-- 有了 star 的函式定義, 可以直接以內定值呼叫 star() 執行繪圖
star()
-- 也可以修改輸入變數, 執行其他繪圖
star(100, 300, 300, 6, "rgb(200,0,0)", "rgb(0,0,200)", "stroke", 45 )

一旦計算機程式進入全球資訊網, 全球協同成為可行模式之一

上述位於 Github 倉儲的 star.lua 程式碼寫法如下:

-- 導入 "js" 模組
local js = require "js"
-- global 就是 javascript 的 window
local global = js.global
local document = global.document
-- html 檔案中 canvas id 設為 "canvas"
local canvas = document:getElementById("canvas")
-- 將 ctx 設為 canvas 2d 繪圖畫布變數
local ctx = canvas:getContext("2d")

-- 乘上 deg 可轉為徑度單位
deg = math.pi / 180

-- 建立多邊形定點位置畫線函式
local function star(radius, xc, yc, n, fs, ss, fors, theta)
    radius = radius or 50
    xc = xc or 100
    yc = yc or 100
    n = n or 5
    -- 屬性呼叫使用 . 而方法呼叫使用 :
    -- 填色屬性
    fs = fs or "rgb(200,0,0)"
    -- 畫筆顏色屬性
    ss = ss or "rgb(0,0,200)"
    -- 內定為填色
    fors = fors or "fill"
    -- 旋轉角度
    theta = theta or 0
    ctx.fillStyle = fs
    ctx.strokeStyle = ss
    xi = xc + radius*math.cos((360/n)*deg+(90+theta)*deg)
    yi = yc - radius*math.sin((360/n)*deg+(90+theta)*deg)
    ctx:beginPath()
    ctx:moveTo(xi,yi)
    for i = 2, n+1 do
        x = xc + radius*math.cos((360/n)*deg*i+(90+theta)*deg)
        y = yc - radius*math.sin((360/n)*deg*i+(90+theta)*deg)
        ctx:lineTo(x,y)
    end
    ctx:closePath()
    if fors == "fill" then
        ctx:fill()
    else
        ctx:stroke()
    end
end

return {
    star = star;
}

概念與知識設計理論說: 概念與知識都是數學理論中的集合

意即, 個人或團隊都可以擁有無限多的概念 (concept) 與無限多的知識 (knowledge), 學習的目的在確認個人或團隊的口袋中 ,到底有哪些概念? 有哪些知識? 簡單說 ,概念就是想法, 而且是當下還沒有知識可以落實的想法, 至於知識, 則是可以用來解決問題, 得到結果的方法集合.

只是, 概念與知識雖然可以無限多, 但工程師解題的過程卻通常有時效限制, 有可用資源的限制, 而且概念與知識都是時間的函式, 也就是說, 先前所獲得的知識, 在空間與時間的遷流變化下, 有可能無法真正落實解題, 而流為概念, 而先前的概念, 在相關知識與實際解題過程, 予以充實後, 可以累積許多有用的知識, 另外, 概念除了透過解題驗證, 可以成為知識外, 也可以觸發出更多概念, 反之, 知識一旦透過證實, 在當下無法解題, 就成為概念之外, 也可以將知識用於解題的過程中, 產生許多其他有用的知識, 相關轉換表示如下:

概念->概念 - 概念可以激發更多概念 知識->知識 - 知識的運用可以獲得更多知識 概念->知識 - 概念的落實當下 ,可以產生知識 知識->概念 - 知識用於解題當下, 若發現已不適用, 則成為概念

段落結論:

假如知識無窮多, 沒有落實不了的概念 產生概念雖不可受限, 但腳踏實地則需要能夠解決問題的知識

可以天馬行空, 同時也要腳踏實地

mdecourse commented 6 years ago

將 Javascript 轉為 web based Lua

<h1>
    yen四連桿</h1>
<script>
var Point = function Point(x,y){
    this.x=x;
    this.y=y;
}
// 附加 drawMe() 方法, 以繪製出點位置
Point.prototype.drawMe = function drawMe(g,r){
    this.g = g;
    this.r = r;
    this.g.save();
    this.g.moveTo(this.x,this.y);
    this.g.beginPath();
    // draw a radius=4 circle
    this.g.arc(this.x, this.y, this.r, 0, 2 * Math.PI, true);
    this.g.moveTo(this.x,this.y);
    this.g.lineTo(this.x+this.r, this.y);
    this.g.moveTo(this.x, this.y);
    this.g.lineTo(this.x-this.r, this.y);
    this.g.moveTo(this.x, this.y);
    this.g.lineTo(this.x, this.y+this.r);
    this.g.moveTo(this.x, this.y);
    this.g.lineTo(this.x, this.y-this.r);
    this.g.restore();
    this.g.stroke();
}
// 加入 Eq 方法
Point.prototype.Eq = function Eq(pt){
    this.x = pt.x;
    this.y = pt.y;
}
// 加入 setPoint 方法
Point.prototype.setPoint = function setPoint(px,py){
    this.x = px;
    this.y = py;
}
// 加上 distance(pt) 方法, 計算點到 pt 的距離
Point.prototype.distance = function distance(pt){
    this.pt = pt;
    return Math.sqrt(Math.pow(this.x-this.pt.x,2)+Math.pow(this.y-this.pt.y,2));
}
// Line 函式物件
var Line = function Line(p1,p2){
    this.p1 = p1;
    this.p2 = p2;
    this.Tail = this.p1;
    this.Head = this.p2;
    this.length = Math.sqrt(Math.pow(this.p2.x-this.p1.x, 2)+Math.pow(this.p2.y-this.p1.y,2));
}
// setPP 方法 for Line
Line.prototype.setPP = function setPP(p1,p2){
    this.p1 = p1;
    this.p2 = p2;
    this.Tail = this.p1;
    this.Head = this.p2;
    this.length = Math.sqrt(Math.pow(this.p2.x-this.p1.x, 2)+Math.pow(this.p2.y-this.p1.y,2));
}
// setRT 方法 for Line, 應該已經確定 Tail 點, 然後以 r, t 作為設定 Head 的參考
Line.prototype.setRT = function setRT(r,t){
    this.r = r;
    this.t = t;
    var x = this.r * Math.cos(this.t);
    var y = this.r * Math.sin(this.t);
    this.Tail.Eq(this.p1);
    this.Head.setPoint(this.Tail.x + x,this.Tail.y + y);
}
// getR 方法 for Line
Line.prototype.getR = function getR(){
    // x 分量與 y 分量
    var x = this.p1.x - this.p2.x;
    var y = this.p1.y - this.p2.y;
    return Math.sqrt(x * x + y * y);
}
// 根據定義 atan2(y,x), 表示 (x,y) 與 正 x 軸之間的夾角, 介於 PI 與 -PI 間
Line.prototype.getT = function getT(){
    var x = this.p2.x - this.p1.x;
    var y = this.p2.y - this.p1.y;
    if (Math.abs(x) < 1E-100) {
        return y < 0.0 ? -Math.PI/2 : Math.PI/2;
    }else{
        return Math.atan2(y, x);
    }
}
// setTail 方法 for Line
Line.prototype.setTail = function setTail(pt){
    this.pt = pt;
    this.Tail.Eq(pt);
    this.Head.setPoint(this.pt.x + this.x, this.pt.y + this.y);
}
// getHead 方法 for Line
Line.prototype.getHead = function getHead(){
    return this.Head;
}
Line.prototype.getTail = function getTail(){
    return this.Tail;
}
Line.prototype.drawMe = function drawMe(g){
    this.g = g;
    this.g.beginPath();
    this.g.moveTo(this.p1.x,this.p1.y);
    this.g.lineTo(this.p2.x,this.p2.y);
    this.g.stroke();
}
// 轉換函式
function degToRad(x) {
    return x / 180 * Math.PI;
}
function radToDeg(x) {
    return x / Math.PI * 180;
}
//
// 建立一個物件繼承函式
// We need a utility function to do the inheritance
function inherit(superClass, subClass) {
    for(var i in superClass.prototype) {
        subClass.prototype[i] = superClass.prototype[i]
    }
}
// 建立 Link function 物件
var Link = function Link(p1,p2){
    this.p1 = p1;
    this.p2 = p2;
    this.length = Math.sqrt(Math.pow(this.p2.x-this.p1.x, 2)+Math.pow(this.p2.y-this.p1.y,2));
}
// 讓 Link 繼承 Line 的方法與屬性
inherit(Line,Link);
//
//
//
// 建立 Link 特有的 drawMe 方法
Link.prototype.drawMe = function drawMe(g){
    this.g = g;
    var hole = 5;
    var radius = 10;
    var length = this.getR();
    // 儲存先前的繪圖狀態
    this.g.save();
    this.g.translate(this.p1.x,this.p1.y);
    // 這裡的轉角必須配合最初的 Link 是畫在 x 軸上或是 y 軸上來進行座標轉換, 目前是以畫在 y 軸上進行座標軸旋轉, 並且確定 Math.atan2(y,x)
    // 以下 alert 用來 debug
    //alert("角度為"+ radToDeg(this.getT()));
    //alert("座標軸轉角為"+radToDeg(-(Math.PI/2-this.getT())));
    this.g.rotate(-(Math.PI/2-this.getT()));
    // 必須配合畫在 y 軸上的 Link, 進行座標轉換, 也可以改為畫在 x 軸上...
    this.g.moveTo(0,0);
    this.g.beginPath();
    this.g.arc(0, 0, hole, 0, 2*Math.PI, true);
    this.g.stroke();
    //
    this.g.moveTo(0,length);
    this.g.beginPath();
    this.g.arc(0,length, hole, 0, 2*Math.PI, true);
    this.g.stroke();
    //
    this.g.moveTo(0,0);
    this.g.beginPath();
    this.g.arc(0,0, radius, 0, Math.PI, true);
    this.g.moveTo(0+radius,0);
    this.g.lineTo(0+radius,0+length);
    this.g.stroke();
    this.g.moveTo(0,0+length);
    //
    this.g.beginPath();
    this.g.arc(0, 0+length, radius, Math.PI, 0, true);
    this.g.moveTo(0-radius,0+length);
    this.g.lineTo(0-radius,0);
    this.g.stroke();
    //
    this.g.restore();
}
//
//
//
//
//
// ap1 角為 p1 點所在的角度, lenp1 長度則為 ap1 角度對應的邊長
// ap2 角為 p2 點所在的角度, lenp2 長度則為 ap2 角度對應的邊長
// ap3 角為 p3 點所在的角度, lenp3 長度則為 ap3 角度對應的邊長
var Triangle = function Triangle(p1,p2,p3){
    // 先將輸入變數轉為函式物件性質
    this.p1 = p1;
    this.p2 = p2;
    this.p3 = p3;
}
//
Triangle.prototype.getLenp3 = function getLenp3(){
    var p1 = this.p1;
    var ret = p1.distance(this.p2);
    return ret;
}
//
Triangle.prototype.getLenp1 = function getLenp1(){
    var p2 = this.p2;
    var ret = p2.distance(this.p3);
    return ret;
}
//
Triangle.prototype.getLenp2 = function getLenp2(){
    var p1 = this.p1;
    var ret = p1.distance(this.p3);
    return ret;
}

// 角度
Triangle.prototype.getAp1 = function getAp1(){
    var ret = Math.acos(((this.getLenp2() * this.getLenp2() + this.getLenp3() * this.getLenp3()) - this.getLenp1() * this.getLenp1()) / (2* this.getLenp2() * this.getLenp3()));
    return ret;
}
//
Triangle.prototype.getAp2 = function getAp2(){
    var ret =Math.acos(((this.getLenp1() * this.getLenp1() + this.getLenp3() * this.getLenp3()) - this.getLenp2() * this.getLenp2()) / (2* this.getLenp1() * this.getLenp3()));
    return ret;
}
//
Triangle.prototype.getAp3 = function getAp3(){
    var ret = Math.acos(((this.getLenp1() * this.getLenp1() + this.getLenp2() * this.getLenp2()) - this.getLenp3() * this.getLenp3()) / (2* this.getLenp1() * this.getLenp2()));
    return ret;
}
//
Triangle.prototype.drawMe = function drawMe(g){
    this.g = g;
    var r = 5;
    // 繪出三個頂點
    this.p1.drawMe(this.g,r);
    this.p2.drawMe(this.g,r);
    this.p3.drawMe(this.g,r);
    var line1 = new Line(this.p1,this.p2);
    var line2 = new Line(this.p1,this.p3);
    var line3 = new Line(this.p2,this.p3);
    // 繪出三邊線
    line1.drawMe(this.g);
    line2.drawMe(this.g);
    line3.drawMe(this.g);
}
// ends Triangle function
// 透過三個邊長定義三角形
Triangle.prototype.setSSS = function setSSS(lenp3,lenp1,lenp2){
    this.lenp3 = lenp3;
    this.lenp1 = lenp1;
    this.lenp2 = lenp2;
    this.ap1 = Math.acos(((this.lenp2 * this.lenp2 + this.lenp3 * this.lenp3) - this.lenp1 * this.lenp1) / (2* this.lenp2 * this.lenp3));
    this.ap2 = Math.acos(((this.lenp1 * this.lenp1 + this.lenp3 * this.lenp3) - this.lenp2 * this.lenp2) / (2* this.lenp1 * this.lenp3));
    this.ap3 = Math.acos(((this.lenp1 * this.lenp1 + this.lenp2 * this.lenp2) - this.lenp3 * this.lenp3) / (2* this.lenp1 * this.lenp2));
}
// ends setSSS
// 透過兩個邊長與夾角定義三角形
Triangle.prototype.setSAS = function setSAS(lenp3,ap2,lenp1){
    this.lenp3 = lenp3;
    this.ap2 = ap2;
    this.lenp1 = lenp1;
    this.lenp2 = Math.sqrt((this.lenp3 * this.lenp3 + this.lenp1 * this.lenp1) - 2* this.lenp3 * this.lenp1 * Math.cos(this.ap2));
    //等於 SSS(AB, BC, CA);
}
// ends setSAS
//
Triangle.prototype.setSaSS = function setSaSS(lenp2,lenp3,lenp1){
    this.lenp2 = lenp2;
    this.lenp3 = lenp3;
    this.lenp1 = lenp1;
    var ret;
    if(this.lenp1 > (this.lenp2 + this.lenp3)){
    // <CAB 夾角為 180 度, 三點共線且 A 介於 BC 之間
        ret = Math.PI;
    } else {
        // <CAB 夾角為 0, 三點共線且 A 不在 BC 之間
        if((this.lenp1 < (this.lenp2 - this.lenp3)) || (this.lenp1 < (this.lenp3 - this.lenp2))){
        ret = 0.0;
        } else {
        // 透過餘絃定理求出夾角 <CAB 
            ret = Math.acos(((this.lenp2 * this.lenp2 + this.lenp3 * this.lenp3) - this.lenp1 * this.lenp1) / (2 * this.lenp2 * this.lenp3));
        }
    }
    return ret;
}
// 取得三角形的三個邊長值
Triangle.prototype.getSSS = function getSSS(){
    var temp = new Array(2);
    temp[0] = this.getLenp1();
    temp[1] = this.getLenp2();
    temp[2] = this.getLenp3();
    return temp;
}
// 取得三角形的三個角度值
Triangle.prototype.getAAA = function getAAA(){
    var temp = new Array(2);
    temp[0] = this.getAp1();
    temp[1] = this.getAp2();
    temp[2] = this.getAp3();
    return temp;
}
// 取得三角形的三個角度與三個邊長
Triangle.prototype.getASASAS = function getASASAS(){
    var temp = new Array(5);
    temp[0] = this.getAp1();
    temp[1] = this.getLenp1();
    temp[2] = this.getAp2();
    temp[3] = this.getLenp2();
    temp[4] = this.getAp3();
    temp[5] = this.getLenp3();
    return temp;
}
Triangle.prototype.setPPSS = function setPPSS(p1,p3,lenp1,lenp3){
var temp = new Array(1);
this.p1 = p1;
this.p3 = p3;
this.lenp1 = lenp1;
this.lenp3 = lenp3;
this.lenp2 = this.p1.distance(this.p3);
// bp3 is the angle beside p3 point, cp3 is the angle for line23, p2 is the output
var ap3,bp3,cp3,p2;
var line31 = new Line(p3,p1);
ap3 = Math.acos(((this.lenp1 * this.lenp1 + this.lenp2 * this.lenp2) - this.lenp3 * this.lenp3) / (2 * this.lenp1 * this.lenp2));
bp3 = line31.getT();
cp3 = bp3 - ap3;
temp[0] = p3.x + this.lenp1*Math.cos(cp3); // p2.x
temp[1] = p3.y + this.lenp1*Math.sin(cp3); // p2.y
return temp;
}
//
//
//
//
// 執行繪圖流程, 注意 x, y 為 global variables
function draw(){
    // 清除畫布
    context.clearRect(0, 0, canvas.width, canvas.height)
    // 畫圖
    line1.drawMe(context);
    line2.drawMe(context);
    line3.drawMe(context);
    // 畫出三角形
    //triangle1.drawMe(context);
    //triangle2.drawMe(context);
    // 旋轉角度進行增量
    theta += dx;
    // 根據旋轉角度計算 p2 點的新位置
    p2.x = p1.x + line1.length*Math.cos(theta*degree);
    p2.y = p1.y - line1.length*Math.sin(theta*degree);
    temp = triangle2.setPPSS(p2,p4,link3_len,link2_len);
    p3.x = temp[0];
    p3.y = temp[1];
}
// ends function draw()
//
//
//
//
// 以上為相關函式物件的定義區
// 全域變數
// 幾何位置輸入變數
var x=10,y=10,r=10;
// 畫布與繪圖內容
var canvas,context;
// 其他輸入變數
var theta = 0;
var degree = Math.PI/180;
var dx = 2;
var dy = 4;
var p1 = new Point(150,100);
var p2 = new Point(150,200);
var p3 = new Point(300,300);
var p4 = new Point(350,100);
var line1 = new Link(p1,p2);
var line2 = new Link(p2,p3);
var line3 = new Link(p3,p4);
var line4 = new Link(p1,p4);
var line5 = new Link(p2,p4);
var link2_len = p2.distance(p3);
var link3_len = p3.distance(p4);
var triangle1 = new Triangle(p1,p2,p4);
var triangle2 = new Triangle(p2,p3,p4);
var temp = new Array(1);
//
//
//
//
// 視窗載入時執行內容
window.onload=function(){
    // 繪圖畫布設定
    canvas = document.getElementById("canvas");
    context = canvas.getContext("2d");
    // 座標轉換, 移動 canvas.height 並且 y 座標變號, 也就是將原點座標移到畫面左下角
    context.translate(0,canvas.height);
    context.scale(1,-1);
    // 座標轉換結束, 之後的繪圖都將在此一座標戲中進行繪製
    /*
    context.translate(dx,dy) <==> context.transform( 1,  0,  0,  1, dx, dy)
    context.rotate(θ)        <==> context.transform( c, -s,  s,  c,  0,  0)
    context.scale(sx,sy)     <==> context.transform(sx,  0,  0, sy,  0,  0)
    */
    //////////
    //////////
    // 以間隔 10 micro seconds 重複呼叫 draw()
    window.setInterval(draw, 10);
}
//
//
//
//
</script>
<p>
    <canvas height="600" id="canvas" style="border:1px solid black;" width="600"></canvas></p>
mdecourse commented 6 years ago

ANSI C Runge-Kutta and gnuplot

/* Runge Kutta for a set of first order differential equations */

#include <stdio.h>
#include <math.h>

#define N 2 /* number of first order equations */
#define dist 0.1 /* stepsize in t*/
#define MAX 30.0 /* max for t */

FILE *output; /* internal filename */
FILE *output1; /* internal filename */
// 利用 pipe 呼叫 gnuplot 繪圖
FILE *pipe;

void runge4(double x, double y[], double step); /* Runge-Kutta function */
double f(double x, double y[], int i); /* function for derivatives */

void main(){

  double t, y[N];
  int j;

  output=fopen("osc.dat", "w"); /* external filename */
  output1=fopen("osc1.dat", "w"); /* external filename */

  y[0]=1.0; /* initial position */
  y[1]=0.0; /* initial velocity */

  //fprintf(output, "0\t%f\n", y[0]);

  for (j=1; j*dist<=MAX ;j++) /* time loop */{

    t=j*dist;
    runge4(t, y, dist);
    fprintf(output, "%f\t%f\n", t, y[0]);
    fprintf(output1, "%f\t%f\n", t, y[1]);
  }

  fclose(output);
  fclose(output1);

  pipe = popen("gnuplot -persist","w");
  //fprintf(pipe,"set term png enhanced font \"v:/fireflysung.ttf\" 18 \n");
  fprintf(pipe,"set term png enhanced font \"y:/wqy-microhei.ttc\" 18 \n");
  //fprintf(pipe,"set yrange [68:70]\n");
  fprintf(pipe,"set output \"test.png\"\n");
  fprintf(pipe, "plot \"osc.dat\" title \"位移\" with lines, \"osc1.dat\" title \"速度\" with lines\n");
  fprintf(pipe,"quit\n");

  fprintf(pipe,"quit\n");
  pclose(pipe);
}

void runge4(double x, double y[], double step){

  double h=step/2.0, /* the midpoint */
  t1[N], t2[N], t3[N], /* temporary storage arrays */
  k1[N], k2[N], k3[N],k4[N]; /* for Runge-Kutta */
  int i;

  for (i=0;i<N;i++){

    t1[i]=y[i]+0.5*(k1[i]=step*f(x,y,i));
  }

  for (i=0;i<N;i++){

    t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
  }

  for (i=0;i<N;i++){

    t3[i]=y[i]+ (k3[i]=step*f(x+h, t2, i));
  }

  for (i=0;i<N;i++){

    k4[i]= step*f(x+step, t3, i);
  }

  for (i=0;i<N;i++){

    y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
  }
}

double f(double x, double y[], int i){

  if (i==0)
    x=y[1]; /* derivative of first equation */
  if (i==1)
    x=-y[0]-0.5*y[1];
  return x;
}

https://sourceforge.net/projects/gnuplot/

wqy-microhei_tcc.zip

mdecourse commented 6 years ago

Lua RK4

-- RungeKutta.lua
-- ----------------------------------------------------------------- --
--      This Lua5 module is Copyright (c) 2010, Peter J Billam       --
--                        www.pjb.com.au                             --
--                                                                   --
--   This module is free software; you can redistribute it and/or    --
--          modify it under the same terms as Lua5 itself.           --
-- ----------------------------------------------------------------- --
local M = {} -- public interface
M.Version = '1.07'
M.VersionDate = '20aug2010'

-- Example usage:
-- local RK = require 'RungeKutta'
-- RK.rk4()

--------------------- infrastructure ----------------------
local function arr2txt(...) -- neat printing of arrays for debug use
    local txt = {}
    for e in ... do txt[#txt+1] = string.format('%g',e) end
    return table.concat(txt,' ') .. "\n"
end
local function warn(str)
    io.stderr:write(str,'\n')
end
local function die(str)
    io.stderr:write(str,'\n')
    os.exit(1)
end
local flag = false
local a
local b
local function gaussn(standdev)
    -- returns normal distribution around 0.0 by the Box-Muller rules
    if not flag then
        a = math.sqrt(-2.0 * math.log(math.random()))
        b = 6.28318531 * math.random()
        flag = true
        return standdev * a * math.sin(b)
    else
        flag = false
        return standdev * a * math.cos(b)
    end
end
-------------------------------------------------------

function M.rk2(yn, dydt, t, dt)
    if type(yn) ~= 'table' then
        warn("RungeKutta.rk2: 1st arg must be an table\n")
        return false
    end
    if type(dydt) ~= 'function' then
        warn("RungeKutta.rk2: 2nd arg must be a function\n")
        return false
    end

    local gamma = .75;  -- Ralston's minimisation of error bounds
    local alpha = 0.5/gamma; local beta = 1.0-gamma;
    local alphadt=alpha*dt; local betadt=beta*dt; local gammadt=gamma*dt;
    local ny = #yn;
    local ynp1 = {}
    local dydtn = {}
    local ynpalpha = {}  -- Gear calls this q
    local dydtnpalpha = {}
    dydtn = dydt(t, yn);
    -- for i=1, ny do
    for i in pairs(yn) do
        ynpalpha[i] = yn[i] + alphadt*dydtn[i];
    end
    dydtnpalpha = dydt(t+alphadt, ynpalpha);
    for i in pairs(yn) do
        ynp1[i] = yn[i]+betadt*dydtn[i]+gammadt*dydtnpalpha[i];
    end
    return t+dt, ynp1
end
function deepcopy(object)  -- http://lua-users.org/wiki/CopyTable
    local lookup_table = {}
    local function _copy(object)
    if type(object) ~= "table" then
        return object
    elseif lookup_table[object] then
        return lookup_table[object]
    end
    local new_table = {}
    lookup_table[object] = new_table
    for index, value in pairs(object) do
        new_table[_copy(index)] = _copy(value)
    end
    return setmetatable(new_table, getmetatable(object))
    end
    return _copy(object)
end

local saved_k0; local use_saved_k0 = false
function M.rk4(yn, dydt, t, dt)
    -- The Runge-Kutta-Merson 5-function-evaluation 4th-order method
    -- in the sine-cosine example, this seems to work as a 7th-order method !
    if (type(yn) ~= 'table') then
        warn("RungeKutta.rk4: 1st arg must be a table\n")
        return false
    end
    if (type(dydt) ~= 'function') then
        warn("RungeKutta.rk4: 2nd arg must be a function\n")
        return false
    end
    local ny = #yn; local i;

    local k0
    if use_saved_k0 then
        k0 = deepcopy(saved_k0)  -- a simpler single-level copy  would do...
-- without the copy() it gets trashed on the 2nd call to this function :-(
    else  k0 = dydt(t, yn)
    end
    for i in pairs(yn) do k0[i] = k0[i] * dt end

    local eta1 = {}
    for i in pairs(yn) do eta1[i] = yn[i] + k0[i]/3.0 end
    local k1 = dydt(t + dt/3.0, eta1)
    for i in pairs(yn) do k1[i] = k1[i] * dt end

    local eta2 = {}
    local k2 = {}
    for i in pairs(yn) do
        eta2[i] = yn[i] + (k0[i]+k1[i])/6.0
    end
    k2 = dydt(t + dt/3.0, eta2)
    for i in pairs(yn) do k2[i] = k2[i] * dt end

    local eta3 = {}
    for i in pairs(yn) do
        eta3[i] = yn[i] + (k0[i]+3.0*k2[i])*0.125
    end
    local k3 = dydt(t+0.5*dt, eta3)
    for i in pairs(yn) do k3[i] = k3[i] * dt end

    local eta4 = {}
    for i in pairs(yn) do
        eta4[i] = yn[i] + (k0[i]-3.0*k2[i]+4.0*k3[i])*0.5
    end
    local k4 = dydt(t+dt, eta4)
    for i in pairs(yn) do k4[i] = k4[i] * dt end

    local ynp1 = {}
    for i in pairs(yn) do
        ynp1[i] = yn[i] + (k0[i]+4.0*k3[i]+k4[i])/6.0;
    end

    -- Merson's method for error estimation, see Gear p85, only works
    -- if F is linear, ie F = Ay + bt, so that eta4 has no 4th-order
    -- errors.  So in general step-doubling is the only way to do it.
    -- Estimate error terms ...
    -- if ($epsilon) {
    --  my $errmax = 0; my $diff;
    --  for ($i=$[; $i<=$ny; $i++) {
    --      $diff = 0.2 * abs ($ynp1[$i] - $eta4[$i]);
    --      if ($errmax < $diff) { $errmax = $diff; }
    --  }
    --  -- print "errmax = $errmax\n"; -- not much related to the actual error
    -- }

    return t+dt, ynp1
end

local t = 0; local halfdt; local y2 = {}
function M.rk4_auto(yn, dydt, t, dt, arg4)
    if (type(yn) ~= 'table') then
        warn("RungeKutta.rk4_auto: 1st arg must be a table\n")
        return false
    end
    if (type(dydt) ~= 'function') then
        warn("RungeKutta.rk4_auto: 2nd arg must be a function\n")
        return false
    end
    if dt == 0 then dt = 0.1 end
    local errors; local epsilon = nil
    if (type(arg4) == 'table') then
        errors = arg4; epsilon = nil
    else
        epsilon = math.abs(arg4); errors = nil
        if epsilon == 0 then epsilon = .0000001 end
    end
    local ny = #yn; local i

    local y1 = {}
    local y3 = {}
    saved_k0 = dydt(t, yn)
    local resizings = 0;
    local highest_low_error = 0.1e-99; local highest_low_dt = 0.0;
    local lowest_high_error = 9.9e99;  local lowest_high_dt = 9.9e99;
    while true do
        halfdt = 0.5 * dt; local dummy
        use_saved_k0 = true
        dummy, y1 = M.rk4(yn, dydt, t, dt)
        dummy, y2 = M.rk4(yn, dydt, t, halfdt)
        use_saved_k0 = false
        dummy, y3 = M.rk4(y2, dydt, t+halfdt, halfdt)

        local relative_error
        if epsilon then
            local errmax = 0; local diff; local ymax = 0
            for i in pairs(yn) do
                diff = math.abs(y1[i] - y3[i])
                if errmax < diff then errmax = diff end
                if ymax < math.abs(yn[i]) then ymax = math.abs(yn[i]) end
            end
            relative_error = errmax / (epsilon*ymax)
        elseif errors then
            relative_error = 0.0; local diff;
            for i in pairs(yn) do
                diff = math.abs(y1[i] - y3[i]) / math.abs(errors[i])
                if relative_error < diff then relative_error = diff end
            end
        else
            die "RungeKutta.rk4_auto: \$epsilon & \@errors both undefined\n";
        end
        -- Gear's "correction" assumes error is always in 5th-order terms :-(
        -- $y1[$i] = (16.0*$y3{$i] - $y1[$i]) / 15.0;
        if relative_error < 0.60 then
            if dt > highest_low_dt then
                highest_low_error = relative_error; highest_low_dt = dt
            end
        elseif relative_error > 1.67 then
            if dt < lowest_high_dt then
                lowest_high_error = relative_error; lowest_high_dt = dt
            end
        else
            break
        end
        if lowest_high_dt<9.8e99 and highest_low_dt>1.0e-99 then -- interpolate
            local denom = math.log(lowest_high_error/highest_low_error)
            if highest_low_dt==0.0 or highest_low_error==0.0 or denom == 0.0 then
                dt = 0.5 * (highest_low_dt+lowest_high_dt)
            else
                dt = highest_low_dt * ( (lowest_high_dt/highest_low_dt)
                 ^ ((math.log(1.0/highest_low_error)) / denom) )
            end
        else
            local adjust = relative_error^(-0.2) -- hope error is 5th-order ...
            if math.abs(adjust) > 2.0 then
                dt = dt * 2.0  -- prevent infinity if 4th-order is exact ...
            else
                dt = dt * adjust
            end
        end
        resizings = resizings + 1
        if resizings>4 and highest_low_dt>1.0e-99 then
            -- hope a small step forward gets us out of this mess ...
            dt = highest_low_dt;  halfdt = 0.5 * dt;
            use_saved_k0 = true
            dummy, y2 = M.rk4(yn, dydt, t, halfdt)
            use_saved_k0 = false
            dummy, y3 = M.rk4(y2, dydt, t+halfdt, halfdt)
            break
        end
    end
    return t+dt, dt, y3
end

function M.rk4_auto_midpoint()
    return t+halfdt, y2
end

------------------------ EXPORT_OK routines ----------------------

function M.rk4_ralston (yn, dydt, t, dt)
    if (type(yn) ~= 'table') then
        warn("RungeKutta.rk4_ralston: 1st arg must be arrayref\n")
        return false
    end
    if (type(dydt) ~= 'function') then
        warn("RungeKutta.rk4_ralston: 2nd arg must be a subroutine ref\n")
        return false
    end
    local ny = #yn; local i;

    -- Ralston's minimisation of error bounds, see Gear p36
    local alpha1=0.4; local alpha2 = 0.4557372542  -- = .875 - .1875*(sqrt 5);

    local k0 = dydt(t, yn)
    for i in pairs(yn) do k0[i] = dt * k0[i] end

    local k1 = {}
    for i in pairs(yn) do k1[i] = yn[i] + 0.4*k0[i] end
    k1 = dydt(t + alpha1*dt, k1)
    for i in pairs(yn) do k1[i] = dt * k1[i] end

    local k2 = {}
    for i in pairs(yn) do
        k2[i] = yn[i] + 0.2969776*k0[i] + 0.15875966*k1[i]
    end
    k2 = dydt(t + alpha2*dt, k2)
    for i in pairs(yn) do k2[i] = dt * k2[i] end

    local k3 = {}
    for i in pairs(yn) do
        k3[i] = yn[i] + 0.21810038*k0[i] - 3.0509647*k1[i] + 3.83286432*k2[i]
    end
    k3 = dydt(t+dt, k3)
    for i in pairs(yn) do k3[i] = dt * k3[i] end

    local ynp1 = {}
    for i in pairs(yn) do
        ynp1[i] = yn[i] + 0.17476028*k0[i]
         - 0.55148053*k1[i] + 1.20553547*k2[i] + 0.17118478*k3[i]
    end
    return t+dt, ynp1
end
function M.rk4_classical(yn, dydt, t, dt)
    if (type(yn) ~= 'table') then
        warn("RungeKutta.rk4_classical: 1st arg must be arrayref\n")
        return false
    end
    if (type(dydt) ~= 'function') then
        warn("RungeKutta.rk4_classical: 2nd arg must be subroutine ref\n")
        return false
    end
    local ny = #yn; local i;

    -- The Classical 4th-order Runge-Kutta Method, see Gear p35

    local k0 = dydt(t, yn)
    for i in pairs(yn) do k0[i] = dt * k0[i] end

    local eta1 = {}
    for i in pairs(yn) do eta1[i] = yn[i] + 0.5*k0[i] end
    local k1 = dydt(t+0.5*dt, eta1)
    for i in pairs(yn) do k1[i] = dt * k1[i] end

    local eta2 = {}
    for i in pairs(yn) do eta2[i] = yn[i] + 0.5*k1[i] end
    local k2 = dydt(t+0.5*dt, eta2)
    for i in pairs(yn) do k2[i] = dt * k2[i] end

    local eta3 = {}
    for i in pairs(yn) do eta3[i] = yn[i] + k2[i] end
    local k3 = dydt(t+dt, eta3)
    for i in pairs(yn) do k3[i] = dt * k3[i] end

    local ynp1 = {}
    for i in pairs(yn) do
        ynp1[i] = yn[i] + (k0[i] + 2.0*k1[i] + 2.0*k2[i] + k3[i]) / 6.0;
    end
    return t+dt, ynp1
end

return M

--[[

=pod

=head1 NAME

RungeKutta.lua - Integrating Systems of Differential Equations

=head1 SYNOPSIS

 local RK = require 'RungeKutta'
 function dydt(t, y) -- the derivative function
   -- y is the table of the values, dydt the table of the derivatives
   -- the table can be an array (1...n), or a dictionary; whichever,
   -- the same indices must be used for the return table: dydt
   local dydt = {}; ... ; return dydt 
 end
 y = initial_y(); t=0; dt=0.4;  -- the initial conditions
 -- For automatic timestep adjustment ...
 while t < tfinal do
    t, dt, y = RK.rk4_auto(y, dydt, t, dt, 0.00001)
    display(t, y)
 end

 -- Or, for fixed timesteps ...
 while t < tfinal do
   t, y = RK.rk4(y, dydt, t, dt)  -- Merson's 4th-order method
   display(t, y)
 end
 -- alternatively, though not so accurate ...
 t, y = RK.rk2(y, dydt, t, dt)   -- Heun's 2nd-order method

 -- or, also available ...
 t, y = RK.rk4_classical(y, dydt, t, dt) -- Runge-Kutta 4th-order
 t, y = RK.rk4_ralston(y, dydt, t, dt)   -- Ralston's 4th-order

=head1 DESCRIPTION

RungeKutta.lua offers algorithms for the numerical integration
of simultaneous differential equations of the form

 dY/dt = F(t,Y)

where Y is an array of variables whose initial values Y(0) are
known, and F is a function known from the dynamics of the problem.

The Runge-Kutta methods all involve evaluating the derivative function
F(t,Y) more than once, at various points within the timestep, and
combining the results to reach an accurate answer for the Y(t+dt).
This module only uses explicit Runge-Kutta methods; the implicit methods
involve, at each timestep, solving a set of simultaneous equations
involving both Y(t) and F(t,Y), and this is generally intractable.

Three main algorithms are offered.  I<rk2> is Heun's 2nd-order
Runge-Kutta algorithm, which is relatively imprecise, but does have
a large range of stability which might be useful in some problems.  I<rk4>
is Merson's 4th-order Runge-Kutta algorithm, which should be the normal
choice in situations where the step-size must be specified.  I<rk4_auto>
uses the step-doubling method to adjust the step-size of I<rk4> automatically
to achieve a specified precision; this saves much fiddling around trying
to choose a good step-size, and can also save CPU time by automatically
increasing the step-size when the solution is changing only slowly.

This module is the translation into I<Lua> of the I<Perl> CPAN module
Math::RungeKutta, and comes in its C<./lua> subdirectory.
There also exists a translation into I<JavaScript>
which comes in its C<./js> subdirectory.
The calling-interfaces are identical in all three versions.

This module has been designed to be robust and easy to use, and should
be helpful in solving systems of differential equations which arise
within a I<Lua> context, such as economic, financial, demographic
or ecological modelling, mechanical or process dynamics, etc.

Version 1.07

=head1 FUNCTIONS

=over 3

=item I<rk2>(y, dydt, t, dt )

where the arguments are:
 I<y> an array of initial values of variables,
 I<dydt> the function calculating the derivatives,
 I<t> the initial time,
 I<dt> the timestep.

The algorithm used is that derived by Ralston, which uses Lotkin's bound
on the derivatives, and minimises the solution error (gamma=3/4).
It is also known as the Heun method, though unfortunately several other
methods are also known under this name. Two function evaluations are needed
per timestep, and the remaining error is in the 3rd and higher order terms.

I<rk2> returns t, y where these are now the new values
at the completion of the timestep.

=item I<rk4>( y, dydt, t, dt )

The arguments are the same as in I<rk2>.

The algorithm used is that developed by Merson,
which performs five function evaluations per timestep.
For the same timestep, I<rk4> is much more accurate than I<rk4_classical>,
so the extra function evaluation is well worthwhile.

I<rk4> returns t, y where these are now the new values
at the completion of the timestep.

=item I<rk4_auto>( y, dydt, t, dt, epsilon )

=item I<rk4_auto>( y, dydt, t, dt, errors )

In the first form the arguments are:
 I<y> an array of initial values of variables,
 I<dydt> the function calculating the derivatives,
 I<t> the initial time,
 I<dt> the initial timestep,
 I<epsilon> the errors per step will be about epsilon*ymax

In the second form the last argument is:
 I<errors> an array of maximum permissible errors.

The first I<epsilon> calling form is useful when all the elements of
I<y> are in the same units and have the same typical size (e.g. y[10]
is population aged 10-11 years, y[25] is population aged 25-26 years).
The default value of the 4th argument is I<epsilon = 0.00001>.

The second I<errors> form is useful otherwise
(e.g. y[1] is gross national product, y[2] is interest rate).
In this calling form, the permissible errors are specified in
absolute size for each variable; they won't get scaled at all.

I<rk4_auto> adjusts the timestep automatically to give the
required precision.  It does this by trying one full-timestep,
then two half-timesteps, and comparing the results.
(Merson's method, as used by I<rk4>, was devised to be able
to give an estimate of the remaining local error; for the
record, it is I<0.2*(ynp1[i]-eta4[i])> in each term.
I<rk4_auto> does not exploit this feature because it only
works for linear I<dydt> functions of the form I<Ay + bt>.)

I<rk4_auto> needs 14 function evaluations per double-timestep, and
it has to re-do 13 of those every time it adjusts the timestep.

I<rk4_auto> returns t, dt, y where these
are now the new values at the completion of the timestep.

=item I<rk4_auto_midpoint>()

I<rk4_auto> performs a double timestep within dt, and returns
the final values; the values as they were at the midpoint do
not normally get returned.  However, if you want to draw a
nice smooth graph, or to update a nice smoothly-moving display,
those values as they were at the midpoint would be useful to you.
Therefore, I<rk4_auto_midpoint> provides a way of retrieving them.

Note that you must call I<rk4_auto> first, which returns the values at
time t+dt, then I<rk4_auto_midpoint> subsequently, which returns the
values at t+dt/2, in other words you get the two sets of values out
of their chronological order. Sorry about this.  For example,

 while t < tfinal do
   t, dt, y = rk4_auto(y, dydt, t, dt, epsilon)
   t_midpoint, y_midpoint = rk4_auto_midpoint()
   update_display(t_midpoint, y_midpoint)
   update_display(t, y)
 end

I<rk4_auto_midpoint> returns t, y where these were the
values at the midpoint of the previous call to I<rk4_auto>.

=back

=head1 CALLER-SUPPLIED FUNCTIONS

=over 3

=item I<dydt>( t, y )

This subroutine will be passed by reference as the second argument to
I<rk2>, I<rk4> and I<rk4_auto>. The name doesn't matter of course.
It must expect the following arguments:
 I<t> the time (in case the equations are time-dependent),
 I<y> the array of values of variables.

It must return an array of the derivatives
of the variables with respect to time.

=back

=head1 EXPORT_OK FUNCTIONS

The following functions are not the usual first choice,
but are supplied in case you need them:

=over 3

=item I<rk4_classical>( y, dydt, t, dt )

The arguments and the return values are the same as in I<rk2> and I<rk4>.

The algorithm used is the classic, elegant, 4th-order Runge-Kutta
method, using four function evaluations per timestep:
 k0 = dt * F(y(n))
 k1 = dt * F(y(n) + 0.5*k0)
 k2 = dt * F(y(n) + 0.5*k1)
 k3 = dt * F(y(n) + k2)
 y(n+1) = y(n) + (k0 + 2*k1 + 2*k2 + k3) / 6

=item I<rk4_ralston>( y, dydt, t, dt )

The arguments and the return values are the same as in I<rk2> and I<rk4>.

The algorithm used is that developed by Ralston, which optimises
I<rk4_classical> to minimise the error bound on each timestep.
This module does not use it as the default 4th-order method I<rk4>,
because Merson's algorithm generates greater accuracy, which allows
the timestep to be increased, which more than compensates for
the extra function evaluation.

=back

=head1 EXAMPLES

There are a couple of example Perl scripts in the I<./examples/>
subdirectory of the build directory.
You can use their code to help you get your first application going.

=over 3

=item I<sine-cosine>

This script uses I<Term::Clui> (arrow keys and Return, or q to quit)
to offer a selection of algorithms, timesteps and error criteria for
the integration of a simple sine/cosine wave around one complete cycle.
This was the script used as a testbed during development.

=item I<three-body>

This script uses the vt100 or xterm 'moveto' and 'reverse'
sequences to display a little simulation of three-body gravity.
It uses I<rk4_auto> because a shorter timestep is needed when
two bodies are close to each other. It also uses I<rk4_auto_midpoint>
to smooth the display.  By changing the initial conditions you
can experience how sensitively the outcome depends on them.

=back

=head1 TRAPS FOR THE UNWARY

Alas, things can go wrong in numerical integration.

One of the most fundamental is B<instability>. If you choose a timestep
I<dt> much larger than time-constants implied in your derivative
function I<dydt>, then the numerical solution will oscillate wildy,
and bear no relation to the real behaviour of the equations.
If this happens, choose a shorter I<dt>.

Some of the most difficult problems involve so-called B<stiff>
derivative functions. These arise when I<dydt> introduces a wide
range of time-constants, from very short to long. In order to avoid
instability, you will have to set I<dt> to correspond to the shortest
time-constant; but this makes it impossibly slow to follow the
evolution of the system over longer times.  You should try to separate
out the long-term part of the problem, by expressing the short-term
process as the finding of some equilibrium, and then assume that that
equilibrium is present and solve the long-term problem on its own.

Similarly, numerical integration doesn't enjoy problems where
time-constants change suddenly, such as balls bouncing off hard
surfaces, etc. You can often tackle these by intervening directly
in the I<@y> array between each timestep. For example, if I<$y[17]>
is the height of the ball above the floor, and I<$y[20]> is the
vertical component of the velocity, do something like

 if y[17]<0.0 then y[17] = -0.9*y[17]; y[20] = -0.9*y[20] end

and thus, again, let the numerical integration solve just the
smooth part of the problem.

=head1 JAVASCRIPT

In the C<js/> subdirectory of the install directory there is I<RungeKutta.js>,
which is an exact translation of this Perl code into JavaScript.
The function names and arguments are unchanged.
Brief Synopsis:

 <SCRIPT type="text/javascript" src="RungeKutta.js"> </SCRIPT>
 <SCRIPT type="text/javascript">
 var dydt = function (t, y) {  // the derivative function
    var dydt_array = new Array(y.length); ... ; return dydt_array;
 }
 var y = new Array();

 // For automatic timestep adjustment ...
 y = initial_y(); var t=0; var dt=0.4;  // the initial conditions
 // Arrays of return vaules:
 var tmp_end = new Array(3);  var tmp_mid = new Array(2);
 while (t < tfinal) {
    tmp_end = rk4_auto(y, dydt, t, dt, 0.00001);
    tmp_mid = rk4_auto_midpoint();
    t=tmp_mid[0]; y=tmp_mid[1];
    display(t, y);   // e.g. could use wz_jsgraphics.js or SVG
    t=tmp_end[0]; dt=tmp_end[1]; y=tmp_end[2];
    display(t, y);
 }

 // Or, for fixed timesteps ...
 y = post_ww2_y(); var t=1945; var dt=1;  // start in 1945
 var tmp = new Array(2);  // Array of return values
 while (t <= 2100) {
    tmp = rk4(y, dydt, t, dt);  // Merson's 4th-order method
    t=tmp[0]; y=tmp[1];
    display(t, y);
 }
 </SCRIPT>

I<RungeKutta.js> uses several global variables
which all begin with the letters C<_rk_> so you should
avoid introducing variables beginning with these characters.

=head1 AUTHOR

Peter J Billam, http://www.pjb.com.au/comp/contact.html

=head1 REFERENCES

I<On the Accuracy of Runge-Kutta's Method>,
M. Lotkin, MTAC, vol 5, pp 128-132, 1951

I<An Operational Method for the study of Integration Processes>,
R. H. Merson,
Proceedings of a Symposium on Data Processing,
Weapons Research Establishment, Salisbury, South Australia, 1957

I<Numerical Solution of Ordinary and Partial Differential Equations>,
L. Fox, Pergamon, 1962

I<A First Course in Numerical Analysis>, A. Ralston, McGraw-Hill, 1965

I<Numerical Initial Value Problems in Ordinary Differential Equations>,
C. William Gear, Prentice-Hall, 1971

=head1 SEE ALSO

See also the scripts examples/sine-cosine and examples/three-body,
http://www.pjb.com.au/,
http://www.pjb.com.au/comp/,
Math::WalshTransform,
Math::Evol,
Term::Clui,
Crypt::Tea_JS,
http://www.xmds.org/

=cut
--]]
mdecourse commented 6 years ago

Windows example

#!/usr/bin/env lua
-- 導入 gd 繪圖程式庫
require "gd"

-- 採用楷書字型
local 字型 = "c:/windows/fonts/kaiu.ttf"

-- 設定繪圖區域大小
繪圖 = gd.createTrueColor(500, 500)

-- 定義繪圖顏色
白色 = 繪圖:colorAllocate(255, 255, 255)
藍色 = 繪圖:colorAllocate(0, 0, 255)

-- 繪圖區域填入白色
繪圖:filledRectangle(0, 0, 500, 500, 白色)

-- 利用環型字串列印方法執行繪圖
-- gdImage:stringFTCircle(cx, cy, radius, textRadius, fillPortion, fontname, points, top, bottom, color)
-- im:stringFTCircle(50, 50, 40, 10, 0.9, "./Vera.ttf", 24, "Powered", "by Lua", blue)
繪圖:stringFTCircle(250, 250, 200, 50, 0.9, 字型, 40, "中文上方文字", "中文下方文字",藍色)

-- 將繪圖區域內容存檔
繪圖:png("out.png")

-- 告知使用者繪圖已經完成
print("已經完成")
-- os.execute("D:\\GIMPPortable\\gimpportable.exe out.png")
mdecourse commented 6 years ago

DE.c

/***************************************************************
**                                                            **
**        D I F F E R E N T I A L     E V O L U T I O N       **
**                                                            **
** Program: de.c                                              **
** Version: 3.6                                               **
**                                                            **
** Authors: Dr. Rainer Storn                                  **
**          c/o ICSI, 1947 Center Street, Suite 600           **
**          Berkeley, CA 94707                                **
**          Tel.:   510-642-4274 (extension 192)              **
**          Fax.:   510-643-7684                              **
**          E-mail: storn@icsi.berkeley.edu                   **
**          WWW: http://http.icsi.berkeley.edu/~storn/        **
**          on leave from                                     **
**          Siemens AG, ZFE T SN 2, Otto-Hahn Ring 6          **
**          D-81739 Muenchen, Germany                         **
**          Tel:    636-40502                                 **
**          Fax:    636-44577                                 **
**          E-mail: rainer.storn@zfe.siemens.de               **
**                                                            **
**          Kenneth Price                                     **
**          836 Owl Circle                                    **
**          Vacaville, CA 95687                               **
**          E-mail: kprice@solano.community.net               ** 
**                                                            **
** This program implements some variants of Differential      **
** Evolution (DE) as described in part in the techreport      **
** tr-95-012.ps of ICSI. You can get this report either via   **
** ftp.icsi.berkeley.edu/pub/techreports/1995/tr-95-012.ps.Z  **
** or via WWW: http://http.icsi.berkeley.edu/~storn/litera.html*
** A more extended version of tr-95-012.ps is submitted for   **
** publication in the Journal Evolutionary Computation.       ** 
**                                                            **
** You may use this program for any purpose, give it to any   **
** person or change it according to your needs as long as you **
** are referring to Rainer Storn and Ken Price as the origi-  **
** nators of the the DE idea.                                 **
** If you have questions concerning DE feel free to contact   **
** us. We also will be happy to know about your experiences   **
** with DE and your suggestions of improvement.               **
**                                                            **
***************************************************************/
/**H*O*C**************************************************************
**                                                                  **
** No.!Version! Date ! Request !    Modification           ! Author **
** ---+-------+------+---------+---------------------------+------- **
**  1 + 3.1  +5/18/95+   -     + strategy DE/rand-to-best/1+  Storn **
**    +      +       +         + included                  +        **
**  1 + 3.2  +6/06/95+C.Fleiner+ change loops into memcpy  +  Storn **
**  2 + 3.2  +6/06/95+   -     + update comments           +  Storn **
**  1 + 3.3  +6/15/95+ K.Price + strategy DE/best/2 incl.  +  Storn **
**  2 + 3.3  +6/16/95+   -     + comments and beautifying  +  Storn **
**  3 + 3.3  +7/13/95+   -     + upper and lower bound for +  Storn **
**    +      +       +         + initialization            +        **
**  1 + 3.4  +2/12/96+   -     + increased printout prec.  +  Storn **
**  1 + 3.5  +5/28/96+   -     + strategies revisited      +  Storn **
**  2 + 3.5  +5/28/96+   -     + strategy DE/rand/2 incl.  +  Storn **
**  1 + 3.6  +8/06/96+ K.Price + Binomial Crossover added  +  Storn **
**  2 + 3.6  +9/30/96+ K.Price + cost variance output      +  Storn **
**  3 + 3.6  +9/30/96+   -     + alternative to ASSIGND    +  Storn **
**  4 + 3.6  +10/1/96+   -    + variable checking inserted +  Storn **
**  5 + 3.6  +10/1/96+   -     + strategy indic. improved  +  Storn **
**                                                                  **
***H*O*C*E***********************************************************/

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "memory.h"
#include <time.h>

// 最大族群數, NP
#define MAXPOP  5000
// 最大向量維度, D
#define MAXDIM  60
// 1 為最大化問題, 0 為最小化問題
#define MAXIMAPROBLEM 0
// 可能要配合最大或最小化進行變號
#define PENALITY 1.0E20

/*------Constants for rnd_uni()--------------------------------------------*/

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

// 與機構合成相關的常數定義
#define PI 3.1415926
#define degree PI/180.0
#define mech_loop -1
#define NUM_OF_POINTS 25

/*------------------------Macros----------------------------------------*/

/*#define ASSIGND(a,b) memcpy((a),(b),sizeof(double)*D) */  /* quick copy by Claudio */
                                                           /* works only for small  */
                                                           /* arrays, but is faster.*/

/*------------------------Globals---------------------------------------*/

long  rnd_uni_init;                 /* serves as a seed for rnd_uni()   */
double c[MAXPOP][MAXDIM], d[MAXPOP][MAXDIM];
double (*pold)[MAXPOP][MAXDIM], (*pnew)[MAXPOP][MAXDIM], (*pswap)[MAXPOP][MAXDIM];

/*---------Function declarations----------------------------------------*/

void  assignd(int D, double a[], double b[]);
double rnd_uni(long *idum);    /* uniform pseudo random number generator */
double extern evaluate(int D, double tmp[], long *nfeval); /* obj. funct. */

// 與機構合成相關的函式宣告
double distance(double x0, double y0, double x1, double y1);
double rr(double L1, double dd, double theta);
struct Coord triangletip_coord(double x0, double y0, double R0, double R1, double x1, double y1, double localt);
void mechanism(double x0, double y0, double x1, double y1, double L1,
  double L2, double L3, double L5, double L6, double input_angles[NUM_OF_POINTS], struct Coord output_points[NUM_OF_POINTS]);
double error_function(struct Coord output_points[NUM_OF_POINTS], struct Coord target_points[NUM_OF_POINTS]);
// 用來利用 tip1 與 tip2 的座標, 以及 r1, r2 求最後的三角形頂點座標
struct Coord finaltip_coord(struct Coord tip1_coord, struct Coord tip2_coord, double r1, double r2);

/*---------Function definitions-----------------------------------------*/
// 指定向量 b 為 a
void  assignd(int D, double a[], double b[])
/**C*F****************************************************************
**                                                                  **
** Assigns D-dimensional vector b to vector a.                      **
** You might encounter problems with the macro ASSIGND on some      **
** machines. If yes, better use this function although it's slower. **
**                                                                  **
***C*F*E*************************************************************/
{
   int j;
   for (j=0; j<D; j++)
   {
      a[j] = b[j];
   }
}

// 產生 0 ~ 1 間的亂數
double rnd_uni(long *idum)
/**C*F****************************************************************
**                                                                  **
** SRC-FUNCTION   :rnd_uni()                                        **
** LONG_NAME      :random_uniform                                   **
** AUTHOR         :(see below)                                      **
**                                                                  **
** DESCRIPTION    :rnd_uni() generates an equally distributed ran-  **
**                 dom number in the interval [0,1]. For further    **
**                 reference see Press, W.H. et alii, Numerical     **
**                 Recipes in C, Cambridge University Press, 1992.  **
**                                                                  **
** FUNCTIONS      :none                                             **
**                                                                  **
** GLOBALS        :none                                             **
**                                                                  **
** PARAMETERS     :*idum    serves as a seed value                  **
**                                                                  **
** PRECONDITIONS  :*idum must be negative on the first call.        **
**                                                                  **
** POSTCONDITIONS :*idum will be changed                            **
**                                                                  **
***C*F*E*************************************************************/
{
  long j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  double temp;

  if (*idum <= 0)
  {
    if (-(*idum) < 1) *idum=1;
    else *idum = -(*idum);
    idum2=(*idum);
    for (j=NTAB+7;j>=0;j--)
    {
      k=(*idum)/IQ1;
      *idum=IA1*(*idum-k*IQ1)-k*IR1;
      if (*idum < 0) *idum += IM1;
      if (j < NTAB) iv[j] = *idum;
    }
    iy=iv[0];
  }
  k=(*idum)/IQ1;
  *idum=IA1*(*idum-k*IQ1)-k*IR1;
  if (*idum < 0) *idum += IM1;
  k=idum2/IQ2;
  idum2=IA2*(idum2-k*IQ2)-k*IR2;
  if (idum2 < 0) idum2 += IM2;
  j=iy/NDIV;
  iy=iv[j]-idum2;
  iv[j] = *idum;
  if (iy < 1) iy += IMM1;
  if ((temp=AM*iy) > RNMX) return RNMX;
  else return temp;

}/*------End of rnd_uni()--------------------------*/

// 將上下限轉為全域變數
double inibound_h;      /* upper parameter bound              */
double inibound_l;      /* lower parameter bound              */
// 與機構合成相關的全域變數
// 宣告一個座標結構
struct Coord {
    double x;
    double y;
  // 這裡保留 double z;
};

int main(int argc, char *argv[])
/**C*F****************************************************************
**                                                                  **
** SRC-FUNCTION   :main()                                           **
** LONG_NAME      :main program                                     **
** AUTHOR         :Rainer Storn, Kenneth Price                      **
**                                                                  **
** DESCRIPTION    :driver program for differential evolution.       **
**                                                                  **
** FUNCTIONS      :rnd_uni(), evaluate(), printf(), fprintf(),      **
**                 fopen(), fclose(), fscanf().                     **
**                                                                  **
** GLOBALS        :rnd_uni_init    input variable for rnd_uni()     **
**                                                                  **
** PARAMETERS     :argc            #arguments = 3                   **
**                 argv            pointer to argument strings      **
**                                                                  **
** PRECONDITIONS  :main must be called with three parameters        **
**                 e.g. like de1 <input-file> <output-file>, if     **
**                 the executable file is called de1.               **
**                 The input file must contain valid inputs accor-  **
**                 ding to the fscanf() section of main().          **
**                                                                  **
** POSTCONDITIONS :main() produces consecutive console outputs and  **
**                 writes the final results in an output file if    **
**                 the program terminates without an error.         **
**                                                                  **
***C*F*E*************************************************************/

{
   char  chr;             /* y/n choice variable                */
   char  *strat[] =       /* strategy-indicator                 */
   {
            "",
            "DE/best/1/exp",
            "DE/rand/1/exp",
            "DE/rand-to-best/1/exp",
            "DE/best/2/exp",
            "DE/rand/2/exp",
            "DE/best/1/bin",
            "DE/rand/1/bin",
            "DE/rand-to-best/1/bin",
            "DE/best/2/bin",
            "DE/rand/2/bin"
   };

   int   i, j, L, n;      /* counting variables                 */
   int   r1, r2, r3, r4;  /* placeholders for random indexes    */
   int   r5;              /* placeholders for random indexes    */
   int   D;               /* Dimension of parameter vector      */
   int   NP;              /* number of population members       */
   int   imin;            /* index to member with lowest energy */
   int   refresh;         /* refresh rate of screen output      */
   int   strategy;        /* choice parameter for screen output */
   int   gen, genmax, seed;   

   long  nfeval;          /* number of function evaluations     */

   double trial_cost;      /* buffer variable   */

   // 將上下限轉為全域變數, 可能要根據各變數加以設定
   //double inibound_h;      /* upper parameter bound              */
   //double inibound_l;      /* lower parameter bound              */
   double tmp[MAXDIM], best[MAXDIM], bestit[MAXDIM]; /* members  */
   double cost[MAXPOP];    /* obj. funct. values                 */
   double cvar;            /* computes the cost variance         */
   double cmean;           /* mean cost                          */
   double F,CR;            /* control variables of DE            */
   double cmin;            /* help variables                     */

   FILE  *fpin_ptr;
   FILE  *fpout_ptr;

// 計算執行過程所需時間起點, 需要導入 time.h
  clock_t start = clock();

/*------Initializations----------------------------*/

 //if (argc != 3)                                 /* number of arguments */
 //{
    //printf("\nUsage : de <input-file> <output-file>\n");
    //exit(1);
 //}

// 將結果寫入 out.dat
 fpout_ptr = fopen("out.dat","w");          /* open output file for reading,    */
                                          /* to see whether it already exists */
 /*
 if ( fpout_ptr != NULL )
 {
    printf("\nOutput file %s does already exist, \ntype y if you ",argv[2]);
    printf("want to overwrite it, \nanything else if you want to exit.\n");
    chr = (char)getchar();
    if ((chr != 'y') && (chr != 'Y'))
    {
      exit(1);
    }
    fclose(fpout_ptr);
 }
*/

/*-----Read input data------------------------------------------------*/

 //fpin_ptr   = fopen(argv[1],"r");
/*
 if (fpin_ptr == NULL)
 {
    printf("\nCannot open input file\n");
    exit(1);
 }*/

 //fscanf(fpin_ptr,"%d",&strategy);       /*---choice of strategy-----------------*/
 //fscanf(fpin_ptr,"%d",&genmax);         /*---maximum number of generations------*/
 //fscanf(fpin_ptr,"%d",&refresh);        /*---output refresh cycle---------------*/
 //fscanf(fpin_ptr,"%d",&D);              /*---number of parameters---------------*/
 //fscanf(fpin_ptr,"%d",&NP);             /*---population size.-------------------*/
 //fscanf(fpin_ptr,"%lf",&inibound_h);    /*---upper parameter bound for init-----*/
 //fscanf(fpin_ptr,"%lf",&inibound_l);    /*---lower parameter bound for init-----*/
 //fscanf(fpin_ptr,"%lf",&F);             /*---weight factor----------------------*/
 //fscanf(fpin_ptr,"%lf",&CR);            /*---crossing over factor---------------*/
 //fscanf(fpin_ptr,"%d",&seed);           /*---random seed------------------------*/
  strategy = 2;
  genmax = 200000;
  // refresh 為每幾筆運算後進行資料列印
  refresh = 100;
  // 配合機構尺寸合成, 每一個體有 9 個機構尺寸值與 25 (NUM_OF_POINTS) 個通過點角度值
  // tmp[0~8] 為機構尺寸, tmp[9~33] 為通過點角度值
  D = 9 + NUM_OF_POINTS;
  NP = 50;
  // 機構變數值上限
  inibound_h = 50.;
  // 機構變數值下限
  inibound_l = 0.;
  // for strategy 1, F=0.9, CR = 1.
  // for strategy 2 F=0.7, CR=0.5
  // 一個小時得到 9.7 的誤差
  // 25 點的題目, 若 penality 只取 1000 則 F = 0.7 似乎為 最大 bound for strategy 1, CR = 1.
  F = 0.85;
  CR = 0.7;
  seed = 3;

 //fclose(fpin_ptr);

/*-----Checking input variables for proper range----------------------------*/

  if (D > MAXDIM)
  {
     printf("\nError! D=%d > MAXDIM=%d\n",D,MAXDIM);
     exit(1);
  }
  if (D <= 0)
  {
     printf("\nError! D=%d, should be > 0\n",D);
     exit(1);
  }
  if (NP > MAXPOP)
  {
     printf("\nError! NP=%d > MAXPOP=%d\n",NP,MAXPOP);
     exit(1);
  }
  if (NP <= 0)
  {
     printf("\nError! NP=%d, should be > 0\n",NP);
     exit(1);
  }
  if ((CR < 0) || (CR > 1.0))
  {
     printf("\nError! CR=%f, should be ex [0,1]\n",CR);
     exit(1);
  }
  if (seed <= 0)
  {
     printf("\nError! seed=%d, should be > 0\n",seed);
     exit(1);
  }
  if (refresh <= 0)
  {
     printf("\nError! refresh=%d, should be > 0\n",refresh);
     exit(1);
  }
  if (genmax <= 0)
  {
     printf("\nError! genmax=%d, should be > 0\n",genmax);
     exit(1);
  }
  if ((strategy < 0) || (strategy > 10))
  {
     printf("\nError! strategy=%d, should be ex {1,2,3,4,5,6,7,8,9,10}\n",strategy);
     exit(1);
  }
  if (inibound_h < inibound_l)
  {
     printf("\nError! inibound_h=%f < inibound_l=%f\n",inibound_h, inibound_l);
     exit(1);
  }

/*-----Open output file-----------------------------------------------*/

   //fpout_ptr   = fopen(argv[2],"w");  /* open output file for writing */

   //if (fpout_ptr == NULL)
   //{
      //printf("\nCannot open output file\n");
      //exit(1);
   //}

/*-----Initialize random number generator-----------------------------*/

 rnd_uni_init = -(long)seed;  /* initialization of rnd_uni() */
 nfeval       =  0;  /* reset number of function evaluations */

/*------Initialization------------------------------------------------*/
/*------Right now this part is kept fairly simple and just generates--*/
/*------random numbers in the range [-initfac, +initfac]. You might---*/
/*------want to extend the init part such that you can initialize-----*/
/*------each parameter separately.------------------------------------*/

   for (i=0; i<NP; i++)
   {
      for (j=0; j<D; j++) /* spread initial population members */
      {
        c[i][j] = inibound_l + rnd_uni(&rnd_uni_init)*(inibound_h - inibound_l);
      }
      cost[i] = evaluate(D,c[i],&nfeval); /* obj. funct. value */
   }
   cmin = cost[0];
   imin = 0;
   for (i=1; i<NP; i++)
   {
     if(MAXIMAPROBLEM == 1)
     {
       // 最大化問題
        if (cost[i]>cmin)
        {
          cmin = cost[i];
          imin = i;
        }
      }
      else
      {
        // 最小化問題
        if (cost[i]<cmin)
        {
          cmin = cost[i];
          imin = i;
        }
      }
   }

   assignd(D,best,c[imin]);            /* save best member ever          */
   assignd(D,bestit,c[imin]);          /* save best member of generation */

   pold = &c; /* old population (generation G)   */
   pnew = &d; /* new population (generation G+1) */

/*=======================================================================*/
/*=========Iteration loop================================================*/
/*=======================================================================*/

   gen = 0;                          /* generation counter reset */
   while ((gen < genmax) /*&& (kbhit() == 0)*/) /* remove comments if conio.h */
   {                                            /* is accepted by compiler    */
      gen++;
      imin = 0;

      for (i=0; i<NP; i++)         /* Start of loop through ensemble  */
      {
     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 2 !!!     */
       r1 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while(r1==i);            

     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 3 !!!     */
       r2 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while((r2==i) || (r2==r1));

     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 4 !!!     */
       r3 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while((r3==i) || (r3==r1) || (r3==r2));

     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 5 !!!     */
       r4 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while((r4==i) || (r4==r1) || (r4==r2) || (r4==r3));

     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 6 !!!     */
       r5 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while((r5==i) || (r5==r1) || (r5==r2) || (r5==r3) || (r5==r4));

/*=======Choice of strategy===============================================================*/
/*=======We have tried to come up with a sensible naming-convention: DE/x/y/z=============*/
/*=======DE :  stands for Differential Evolution==========================================*/
/*=======x  :  a string which denotes the vector to be perturbed==========================*/
/*=======y  :  number of difference vectors taken for perturbation of x===================*/
/*=======z  :  crossover method (exp = exponential, bin = binomial)=======================*/
/*                                                                                        */
/*=======There are some simple rules which are worth following:===========================*/
/*=======1)  F is usually between 0.5 and 1 (in rare cases > 1)===========================*/
/*=======2)  CR is between 0 and 1 with 0., 0.3, 0.7 and 1. being worth to be tried first=*/
/*=======3)  To start off NP = 10*D is a reasonable choice. Increase NP if misconvergence=*/
/*           happens.                                                                     */
/*=======4)  If you increase NP, F usually has to be decreased============================*/
/*=======5)  When the DE/best... schemes fail DE/rand... usually works and vice versa=====*/

/*=======EXPONENTIAL CROSSOVER============================================================*/

/*-------DE/best/1/exp--------------------------------------------------------------------*/
/*-------Our oldest strategy but still not bad. However, we have found several------------*/
/*-------optimization problems where misconvergence occurs.-------------------------------*/
// 1 為最原始的解題邏輯方法
     if (strategy == 1) /* strategy DE0 (not in our paper) */
     {
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D);
       L = 0;
       do
       {                       
         tmp[n] = bestit[n] + F*((*pold)[r2][n]-(*pold)[r3][n]);
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }
/*-------DE/rand/1/exp-------------------------------------------------------------------*/
/*-------This is one of my favourite strategies. It works especially well when the-------*/
/*-------"bestit[]"-schemes experience misconvergence. Try e.g. F=0.7 and CR=0.5---------*/
/*-------as a first guess.---------------------------------------------------------------*/
  // 配合邏輯方法 2 選用 R=0.7, CR=0.5
   else if (strategy == 2) /* strategy DE1 in the techreport */
     {
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D);
       L = 0;
       do
       {                       
         tmp[n] = (*pold)[r1][n] + F*((*pold)[r2][n]-(*pold)[r3][n]);
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }
/*-------DE/rand-to-best/1/exp-----------------------------------------------------------*/
/*-------This strategy seems to be one of the best strategies. Try F=0.85 and CR=1.------*/
/*-------If you get misconvergence try to increase NP. If this doesn't help you----------*/
/*-------should play around with all three control variables.----------------------------*/
     // 方法 3 建議 F=0.85 CR=1.
   else if (strategy == 3) /* similiar to DE2 but generally better */
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
       L = 0;
       do
       {                       
         tmp[n] = tmp[n] + F*(bestit[n] - tmp[n]) + F*((*pold)[r1][n]-(*pold)[r2][n]);
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }
/*-------DE/best/2/exp is another powerful strategy worth trying--------------------------*/
     else if (strategy == 4)
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
       L = 0;
       do
       {                           
         tmp[n] = bestit[n] + 
              ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F;
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }
/*-------DE/rand/2/exp seems to be a robust optimizer for many functions-------------------*/
     else if (strategy == 5)
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
       L = 0;
       do
       {                           
         tmp[n] = (*pold)[r5][n] + 
              ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F;
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }

/*=======Essentially same strategies but BINOMIAL CROSSOVER===============================*/

/*-------DE/best/1/bin--------------------------------------------------------------------*/
     else if (strategy == 6) 
     {
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = bestit[n] + F*((*pold)[r2][n]-(*pold)[r3][n]);
         }
         n = (n+1)%D;
           }
     }
/*-------DE/rand/1/bin-------------------------------------------------------------------*/
     else if (strategy == 7) 
     {
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = (*pold)[r1][n] + F*((*pold)[r2][n]-(*pold)[r3][n]);
         }
         n = (n+1)%D;
           }
     }
/*-------DE/rand-to-best/1/bin-----------------------------------------------------------*/
     else if (strategy == 8) 
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = tmp[n] + F*(bestit[n] - tmp[n]) + F*((*pold)[r1][n]-(*pold)[r2][n]);
         }
         n = (n+1)%D;
           }
     }
/*-------DE/best/2/bin--------------------------------------------------------------------*/
     else if (strategy == 9)
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = bestit[n] + 
              ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F;
         }
         n = (n+1)%D;
           }
     }
/*-------DE/rand/2/bin--------------------------------------------------------------------*/
     else
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = (*pold)[r5][n] + 
              ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F;
         }
         n = (n+1)%D;
           }
     }

/*=======Trial mutation now in tmp[]. Test how good this choice really was.==================*/
     trial_cost = evaluate(D,tmp,&nfeval);  /* Evaluate new vector in tmp[] */
   if(MAXIMAPROBLEM == 1)
   {
    // 改為最大化
       if (trial_cost >= cost[i])   /* improved objective function value ? */
       {                                  
          cost[i]=trial_cost;         
          assignd(D,(*pnew)[i],tmp);
          if (trial_cost>cmin)          /* Was this a new minimum? */
          {                               /* if so...*/
             cmin=trial_cost;           /* reset cmin to new low...*/
             imin=i;
             assignd(D,best,tmp);           
          }                           
       }                            
       else
       {
          assignd(D,(*pnew)[i],(*pold)[i]); /* replace target with old value */
       }
    }
    else
    {
          // 最小化問題
       if (trial_cost <= cost[i])   /* improved objective function value ? */
       {                                  
          cost[i]=trial_cost;         
          assignd(D,(*pnew)[i],tmp);
          if (trial_cost<cmin)          /* Was this a new minimum? */
          {                               /* if so...*/
             cmin=trial_cost;           /* reset cmin to new low...*/
             imin=i;
             assignd(D,best,tmp);           
          }                           
       }                            
       else
       {
          assignd(D,(*pnew)[i],(*pold)[i]); /* replace target with old value */
       }
    }

      }   /* End mutation loop through pop. */

      assignd(D,bestit,best);  /* Save best population member of current iteration */

      /* swap population arrays. New generation becomes old one */

      pswap = pold;
      pold  = pnew;
      pnew  = pswap;

/*----Compute the energy variance (just for monitoring purposes)-----------*/

      cmean = 0.;          /* compute the mean value first */
      for (j=0; j<NP; j++)
      {
         cmean += cost[j];
      }
      cmean = cmean/NP;

      cvar = 0.;           /* now the variance              */
      for (j=0; j<NP; j++)
      {
         cvar += (cost[j] - cmean)*(cost[j] - cmean);
      }
      cvar = cvar/(NP-1);

/*----Output part----------------------------------------------------------*/

      if (gen%refresh==1)   /* display after every refresh generations */
      { /* ABORT works only if conio.h is accepted by your compiler */
    printf("\n\n                         PRESS ANY KEY TO ABORT"); 
    printf("\n\n\n Best-so-far cost funct. value=%-15.10g\n",cmin);

    for (j=0;j<D;j++)
    {
      printf("\n best[%d]=%-15.10g",j,best[j]);
    }
    printf("\n\n Generation=%d  NFEs=%ld   Strategy: %s    ",gen,nfeval,strat[strategy]);
    printf("\n NP=%d    F=%-4.2g    CR=%-4.2g   cost-variance=%-10.5g\n",
               NP,F,CR,cvar);
      }

      fprintf(fpout_ptr,"%ld   %-15.10g\n",nfeval,cmin);
   }
/*=======================================================================*/
/*=========End of iteration loop=========================================*/
/*=======================================================================*/

/*-------Final output in file-------------------------------------------*/

   fprintf(fpout_ptr,"\n\n\n Best-so-far obj. funct. value = %-15.10g\n",cmin);

   for (j=0;j<D;j++)
   {
     fprintf(fpout_ptr,"\n best[%d]=%-15.10g",j,best[j]);
   }
   fprintf(fpout_ptr,"\n\n Generation=%d  NFEs=%ld   Strategy: %s    ",gen,nfeval,strat[strategy]);
   fprintf(fpout_ptr,"\n NP=%d    F=%-4.2g    CR=%-4.2g    cost-variance=%-10.5g\n",
           NP,F,CR,cvar); 

  fclose(fpout_ptr);

  /* Code you want timed here */
  printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
  return(0);
}

/*-----------End of main()------------------------------------------*/

// 適應函式 fittness function (cost function)
double evaluate(int D, double tmp[], long *nfeval)
{
  // 先處理通過 5 個點的四連桿問題
  // x0, y0 為左方固定點座標, 必須在 0 ~ 100 間 - 設為 tmp[0], tmp[1]
  // x1, y1 為右方固定點座標, 必須在 0 ~ 100 間 - 設為 tmp[2], tmp[3]
  // L1 為第一桿件的長度, 必須 > 0, 且小於 100 - 設為 tmp[4]
  // L2 為第二桿件的長度, 必須 > 0, 且小於 100 - 設為 tmp[5]
  // L3 為第三桿件的長度, 必須 > 0, 且小於 100 - 設為 tmp[6]
  // L5, L6 與 L2 共同圍出可動桿 L2 對應的三角形, 關注的三角形頂點即 L5 與 L6 的交點, 而 angle3 則為 L6 之對應角(為固定值)
  // L5, L6 必須 > 0, 且小於 100 - 設為 tmp[7], tmp[8]
  // 以下的角度輸入值, 會隨著目標點數的增加而增加, 其索引值由 9 + 通過點數 - 1 決定, 5 點, 則索引至 13, 若通過 25 點, 則索引值為 9 + 24 = 33
  // input_angles[] 為五的輸入的雙浮點數角度值,代表個體的角度向量值 - 分別設為 tmp[9], tmp[10], tmp[11], tmp[12], tmp[13]
  // 這裡的輸入角度值, 將採用以第一角度>0 作為起點, 隨後則為角度增量, 也都必須大於零
  // output_points[] 則為與 input_angles[] 對應的五個三角形的頂點值, 為座標結構值, 分別有 x 與 y 分量值
  // 當利用個體的向量值, 代入 mechanism 後所得到得 output_points[] 再與 target_points[] 進行 cost function 的誤差值最小化
  /* void mechanism(double x0, double y0, double x1, double y1, double L1,
  double L2, double L3, double L5, double L6, double input_angles[NUM_OF_POINTS], struct Coord output_points[NUM_OF_POINTS])*/
  struct Coord target_points[NUM_OF_POINTS], output_points[NUM_OF_POINTS];
  double input_angles[NUM_OF_POINTS], result, total_angle;
  int i;

  (*nfeval)++;

  target_points[0].x = 4.5;
  target_points[0].y = 6.75;

  target_points[1].x = 5.07;
  target_points[1].y = 6.85;

  target_points[2].x = 5.45;
  target_points[2].y = 6.84;

  target_points[3].x = 5.89;
  target_points[3].y = 6.83;

  target_points[4].x = 6.41;
  target_points[4].y = 6.8;

  target_points[5].x = 6.92;
  target_points[5].y = 6.58;

  target_points[6].x = 7.03;
  target_points[6].y = 5.99;

  target_points[7].x = 6.95;
  target_points[7].y = 5.45;

  target_points[8].x = 6.77;
  target_points[8].y = 5.03;

  target_points[9].x = 6.4;
  target_points[9].y = 4.6;

  target_points[10].x = 5.91;
  target_points[10].y = 4.03;

  target_points[11].x = 5.43;
  target_points[11].y = 3.56;

  target_points[12].x = 4.93;
  target_points[12].y = 2.94;

  target_points[13].x = 4.67;
  target_points[13].y = 2.6;

  target_points[14].x = 4.38;
  target_points[14].y = 2.2;

  target_points[15].x = 4.04;
  target_points[15].y = 1.67;

  target_points[16].x = 3.76;
  target_points[16].y = 1.22;

  target_points[17].x = 3.76;
  target_points[17].y = 1.97;

  target_points[18].x = 3.76;
  target_points[18].y = 2.78;

  target_points[19].x = 3.76;
  target_points[19].y = 3.56;

  target_points[20].x = 3.76;
  target_points[20].y = 4.34;

  target_points[21].x = 3.76;
  target_points[21].y = 4.91;

  target_points[22].x = 3.76;
  target_points[22].y = 5.47;

  target_points[23].x = 3.8;
  target_points[23].y = 5.98;

  target_points[24].x = 4.07;
  target_points[24].y = 6.4;

  // 輸入角度值與 tmp[] 的設定
  for(i = 0; i < NUM_OF_POINTS; i++)
  {
    input_angles[i] = tmp[i + 9];
  }

  // 呼叫 mechanism() 以便計算 output_points[]
  mechanism(tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], tmp[6], tmp[7], tmp[8], input_angles, output_points);

  // for debug
  /*
  if(*nfeval%5000 == 0)
  {
    for(i = 0; i < NUM_OF_POINTS; i++)
    {
      printf("%-15.10g : %-15.10g\n", output_points[i].x, output_points[i].y);
    }
    printf("#####################################\n");
  }
*/
  // 利用 output_points[] 與 target_points 計算誤差值, 該誤差值就是 cost
  result = error_function(output_points, target_points);
  // 這裡要分別針對各變數的約束條件與範圍值來設定傳回 PENALITY 或 誤差值 result

 // x0 與 x1 點位於 -50 與 50 中間
    for(i = 0; i < 4; i++)
  {
    if(tmp[i] < -100 || tmp[i] > 100){
      return PENALITY;
    }
  }

  // 三個連桿值, 一定要為正
    for(i = 4; i < 7; i++)
  {
    if(tmp[i] < 0 || tmp[i] > 100){
      return PENALITY;
    }
  }

    // L5 L6 可以為 0 或負值
    for(i = 7; i < 9; i++)
  {
    if(tmp[i] < -100 || tmp[i] > 100){
      return PENALITY;
    }
  }

  // 角度值不可以小於 0

  for(i = 1; i <= NUM_OF_POINTS; i++)
  {
    if(tmp[D-i] < 0){
      return PENALITY;
    }
  }

  /*
  for(i = 0; i < D - NUM_OF_POINTS; i++)
  {
      if((tmp[i] <= inibound_l)|| (tmp[i] >inibound_h)){
      return PENALITY;
    }
  }

  for(i = 1; i <= NUM_OF_POINTS; i++)
  {
    if(tmp[D-i] < 0){
      return PENALITY;
    }
  }
  */
  return result;

  /*
   double result=0, surface = 80.0, z, volume, penality;
   (*nfeval)++;
   z = (surface-tmp[0]*tmp[1])/(2.0*(tmp[0]+tmp[1]));
   volume = tmp[0]*tmp[1]*z;

  if(volume <= 0){
    return PENALITY;
  }

  if((tmp[0] <= inibound_l)|| (tmp[0] >inibound_h)){
    return PENALITY;
  }

  if((tmp[1] <= inibound_l) || (tmp[1] >inibound_h)){
    return PENALITY;
  }
  // volume must >0 and max volume
  // 目前為最小化問題
   return 1+1/(volume*volume);
   */
}

struct Coord triangletip_coord( double x0, double y0, double R0, double R1, double x1, double y1, double localt)
{
    struct Coord tip_coord;

    if (localt>=0 && localt <PI)
    {
        // 目前蓋掉的式子為利用手動代換出來的版本
        //x_value = ((x1-x0)*(-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2))/(2*sqrt(pow((y1-y0),2)+pow((x1-x0),2)))-(y1-y0)*(-mech_loop)*sqrt(fabs(pow(R0,2)-pow((-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2)),2)/(4*(pow((y1-y0),2)+pow((x1-x0),2))))))/sqrt(pow((y1-y0),2)+pow((x1-x0),2))+x0;
        // 以下的式子,先利用文字編輯器,將原先 stringout() 出來的  sqrt 替換成  sqrtt, 以防止被  maxima 中的 subst("^"=pow,expr) 所替換, subst 之後,再使用文字編輯器換回來,就可以得到正確的 C 對應碼.
        tip_coord.x = pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*(mech_loop*(y1-y0)*sqrt(fabs(pow(R0,2)-pow(pow(y1-y0,2)+
    pow(x1-x0,2),-1)*pow(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2),2)/4))+(x1-x0)*pow(sqrt(pow(y1-y0,2)+
    pow(x1-x0,2)),-1)*(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2))/2)+x0;
    }
    else
    {
        // 目前蓋掉的式子為利用手動代換出來的版本
        //x_value = ((x1-x0)*(-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2))/(2*sqrt(pow((y1-y0),2)+pow((x1-x0),2)))-(y1-y0)*(mech_loop)*sqrt(fabs(pow(R0,2)-pow((-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2)),2)/(4*(pow((y1-y0),2)+pow((x1-x0),2))))))/sqrt(pow((y1-y0),2)+pow((x1-x0),2))+x0;
        tip_coord.x = pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*(-mech_loop*(y1-y0)*sqrt(fabs(pow(R0,2)-pow(pow(y1-y0,2)+
    pow(x1-x0,2),-1)*pow(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2),2)/4))+(x1-x0)*pow(sqrt(pow(y1-y0,2)+
    pow(x1-x0,2)),-1)*(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+
    pow(x1-x0,2))/2)+x0;
    }

// 請注意,與 Maxma 公式中的差異為,在 sqrt()中加入 fabs(),避免因為sqrt()中的負值而造成 NaN (Not a Number 問題.
    if (localt>=0 && localt <PI)
    {
        tip_coord.y = /*((x1-x0)*(-mech_loop)*sqrt(
                fabs(pow(R0,2)-pow((-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2)),2)
                /(4*(pow((y1-y0),2)+pow((x1-x0),2)))
                ))
                +(y1-y0)*(-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2))/(2*sqrt(pow((y1-y0),2)+pow((x1-x0),2))))/sqrt(pow((y1-y0),2)+pow((x1-x0),2))
                +y0;*/
                // 利用 sqrtt 居中進行代換所得到的式子
                pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*((y1-y0)*pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*(-pow(R1,2)+
    pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2))/2-mech_loop*(x1-x0)*sqrt(fabs(pow(R0,2)-pow(pow(y1-y0,2)+
    pow(x1-x0,2),-1)*pow(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2),2)/4)))+y0;

    }
    else
    {
        tip_coord.y = /*((x1-x0)*(mech_loop)*sqrt(
                fabs(pow(R0,2)-pow((-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2)),2)
                /(4*(pow((y1-y0),2)+pow((x1-x0),2)))
                ))
                +(y1-y0)*(-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2))/(2*sqrt(pow((y1-y0),2)+pow((x1-x0),2))))/sqrt(pow((y1-y0),2)+pow((x1-x0),2))
                +y0;*/
                pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*((y1-y0)*pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*(-pow(R1,2)+
    pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2))/2+mech_loop*(x1-x0)*sqrt(fabs(pow(R0,2)-pow(pow(y1-y0,2)+
    pow(x1-x0,2),-1)*pow(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2),2)/4)))+y0;
    }

  return tip_coord;
}

double distance(double x0, double y0, double x1, double y1)
{
    double distance_value;
    distance_value = sqrt(pow((x1-x0),2) + pow((y1-y0),2));
    return distance_value;
}

double rr(double L1, double dd, double theta)
{
    double rr_value;
    rr_value = sqrt(L1*L1+dd*dd - 2*L1*dd*cos(theta));
    return rr_value;
}

// 輸入每一個體的變數向量, 然後求各三角形頂點的座標陣列[NUM_OF_POINTS]
void mechanism(double x0, double y0, double x1, double y1, double L1,
  double L2, double L3, double L5, double L6, double input_angles[NUM_OF_POINTS], struct Coord output_points[NUM_OF_POINTS])
{
  // 這裡的輸入角度, 改採第一角度為起始角, 隨後的角度值則為增量值
  // 此函式要輸入控制變數, 然後計算機構尺寸合成的關注點座標位置
  // 以下為可能的處理變數宣告
  // 這裡希望能夠定義一個 struct 來處理座標點
  double rr_length, dd_length, angle;
  struct Coord link1_tip, link2_tip, triangle_tip;
    double angle2, angle3;
  int i;

  // 開始進行三角形頂點座標的計算
  // 以下變數由每一個體向量提供
  /*
    x0 = 0.0;
    y0 = 0.0;
    x1 = 10.0;
    y1 = 0.0;
    L1 = 5.0;
    L2 = 20;
    L3 = 10;
    L5 = 10;
    L6 = 10;
  */
  dd_length = distance(x0, y0, x1, y1);
  /* 設法表示 triangle 所對應的 local 角度,表示為已知變數與 t 的函式 */
  // 假如採用 finaltip_coord, 則不需要 angle3
  angle3 = acos((pow(L2,2)+pow(L5,2)-pow(L6,2))/(2*L2*L5));

  for(i = 0; i < NUM_OF_POINTS; i++)
  {
    // 先建立第一點座標, 即 i=0 者
    // i=0;
    // angle = i*degree;

    /*
    // 利用角度增量進行運算, 相對於 input_angles[0] 作為基準
    if(i > 0)
    {
      input_angles[i] = input_angles[i] + input_angles[i-1];
    }
    */
    angle = input_angles[i]*degree;
    rr_length = rr(L1, dd_length, angle);
    // 第一次三角形疊代
    link1_tip = triangletip_coord(x0, y0, L1, rr_length, x1, y1, angle);
    // 第二次三角形疊代
    /* 設法表示 link2 所對應的 local 角度,表示為已知變數與 t 的函式 */
    angle2 = acos((pow(L2,2)+pow(rr_length,2)-pow(L3,2))/(2*L2*rr_length));
    link2_tip = triangletip_coord(link1_tip.x, link1_tip.y, L2, L3, x1, y1, angle2);
    // 第三次三角形疊代
    //triangle_tip = triangletip_coord(link1_tip.x, link1_tip.y, L5, L6, link2_tip.x, link2_tip.y, angle3);
    // output_points[i] = triangletip_coord(link1_tip.x, link1_tip.y, L5, L6, link2_tip.x, link2_tip.y, angle3);
    // 這裡要嘗試利用 finaltip_coord() 求 tip3 座標, 而 L5 與 L6 可 0 可負
    output_points[i] = finaltip_coord(link1_tip, link2_tip, L5, L6);
  }
}

double error_function(struct Coord output_points[NUM_OF_POINTS], struct Coord target_points[NUM_OF_POINTS])
{
  double error = 0.0;
  int i;
  for(i = 0; i < NUM_OF_POINTS; i++)
  {
    error += fabs(distance(output_points[i].x, output_points[i].y, target_points[i].x, target_points[i].y));
  }
  return error;
}

struct Coord finaltip_coord(struct Coord tip1_coord, struct Coord tip2_coord, double r1, double r2)
{
  struct Coord tip3_coord;
  double theta3, theta4, length3, length4;
  length3 = sqrt(pow(tip2_coord.x - tip1_coord.x,2) + pow(tip2_coord.y - tip1_coord.y,2));
  length4 = sqrt(pow(r1,2) + pow(r2,2));  
  theta3 = acos((tip2_coord.x - tip1_coord.x) / length3);
  theta4 = acos(r1 / length4);
  tip3_coord.x = tip1_coord.x + length4 * cos(theta3 + theta4);
  tip3_coord.y = tip1_coord.y + length4 * sin(theta3 + theta4);

  return tip3_coord;
}
mdecourse commented 6 years ago

DE10.c

// 必須在演算過程中, 設法限制各變數的上下限!!! 否則演化非常容易發散??

/***************************************************************
**                                                            **
**        D I F F E R E N T I A L     E V O L U T I O N       **
**                                                            **
** Program: de.c                                              **
** Version: 3.6                                               **
**                                                            **
** Authors: Dr. Rainer Storn                                  **
**          c/o ICSI, 1947 Center Street, Suite 600           **
**          Berkeley, CA 94707                                **
**          Tel.:   510-642-4274 (extension 192)              **
**          Fax.:   510-643-7684                              **
**          E-mail: storn@icsi.berkeley.edu                   **
**          WWW: http://http.icsi.berkeley.edu/~storn/        **
**          on leave from                                     **
**          Siemens AG, ZFE T SN 2, Otto-Hahn Ring 6          **
**          D-81739 Muenchen, Germany                         **
**          Tel:    636-40502                                 **
**          Fax:    636-44577                                 **
**          E-mail: rainer.storn@zfe.siemens.de               **
**                                                            **
**          Kenneth Price                                     **
**          836 Owl Circle                                    **
**          Vacaville, CA 95687                               **
**          E-mail: kprice@solano.community.net               ** 
**                                                            **
** This program implements some variants of Differential      **
** Evolution (DE) as described in part in the techreport      **
** tr-95-012.ps of ICSI. You can get this report either via   **
** ftp.icsi.berkeley.edu/pub/techreports/1995/tr-95-012.ps.Z  **
** or via WWW: http://http.icsi.berkeley.edu/~storn/litera.html*
** A more extended version of tr-95-012.ps is submitted for   **
** publication in the Journal Evolutionary Computation.       ** 
**                                                            **
** You may use this program for any purpose, give it to any   **
** person or change it according to your needs as long as you **
** are referring to Rainer Storn and Ken Price as the origi-  **
** nators of the the DE idea.                                 **
** If you have questions concerning DE feel free to contact   **
** us. We also will be happy to know about your experiences   **
** with DE and your suggestions of improvement.               **
**                                                            **
***************************************************************/
/**H*O*C**************************************************************
**                                                                  **
** No.!Version! Date ! Request !    Modification           ! Author **
** ---+-------+------+---------+---------------------------+------- **
**  1 + 3.1  +5/18/95+   -     + strategy DE/rand-to-best/1+  Storn **
**    +      +       +         + included                  +        **
**  1 + 3.2  +6/06/95+C.Fleiner+ change loops into memcpy  +  Storn **
**  2 + 3.2  +6/06/95+   -     + update comments           +  Storn **
**  1 + 3.3  +6/15/95+ K.Price + strategy DE/best/2 incl.  +  Storn **
**  2 + 3.3  +6/16/95+   -     + comments and beautifying  +  Storn **
**  3 + 3.3  +7/13/95+   -     + upper and lower bound for +  Storn **
**    +      +       +         + initialization            +        **
**  1 + 3.4  +2/12/96+   -     + increased printout prec.  +  Storn **
**  1 + 3.5  +5/28/96+   -     + strategies revisited      +  Storn **
**  2 + 3.5  +5/28/96+   -     + strategy DE/rand/2 incl.  +  Storn **
**  1 + 3.6  +8/06/96+ K.Price + Binomial Crossover added  +  Storn **
**  2 + 3.6  +9/30/96+ K.Price + cost variance output      +  Storn **
**  3 + 3.6  +9/30/96+   -     + alternative to ASSIGND    +  Storn **
**  4 + 3.6  +10/1/96+   -    + variable checking inserted +  Storn **
**  5 + 3.6  +10/1/96+   -     + strategy indic. improved  +  Storn **
**                                                                  **
***H*O*C*E***********************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <time.h>

// 最大族群數, NP
#define MAXPOP  5000
// 最大向量維度, D
#define MAXDIM  35
#define MAXIMAPROBLEM 0
#define PENALITY 1000

/*------Constants for rnd_uni()--------------------------------------------*/

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

// 與機構合成相關的常數定義
#define PI 3.1415926
#define degree PI/180.0
#define mech_loop -1
#define NUM_OF_POINTS 10

/*------------------------Macros----------------------------------------*/

/*#define ASSIGND(a,b) memcpy((a),(b),sizeof(double)*D) */  /* quick copy by Claudio */
                                                           /* works only for small  */
                                                           /* arrays, but is faster.*/

/*------------------------Globals---------------------------------------*/

long  rnd_uni_init;                 /* serves as a seed for rnd_uni()   */
double c[MAXPOP][MAXDIM], d[MAXPOP][MAXDIM];
double (*pold)[MAXPOP][MAXDIM], (*pnew)[MAXPOP][MAXDIM], (*pswap)[MAXPOP][MAXDIM];

/*---------Function declarations----------------------------------------*/

void  assignd(int D, double a[], double b[]);
double rnd_uni(long *idum);    /* uniform pseudo random number generator */
double extern evaluate(int D, double tmp[], long *nfeval); /* obj. funct. */

// 與機構合成相關的函式宣告
double distance(double x0, double y0, double x1, double y1);
double rr(double L1, double dd, double theta);
struct Coord triangletip_coord( double x0, double y0, double R0, double R1, double x1, double y1, double localt);
void mechanism(double x0, double y0, double x1, double y1, double L1,
  double L2, double L3, double L5, double L6, double input_angles[NUM_OF_POINTS], struct Coord output_points[NUM_OF_POINTS]);
double error_function(struct Coord output_points[NUM_OF_POINTS], struct Coord target_points[NUM_OF_POINTS]);

struct Coord finaltip_coord(struct Coord tip1_coord, struct Coord tip2_coord, double r1, double r2);

/*---------Function definitions-----------------------------------------*/
// 指定向量 b 為 a
void  assignd(int D, double a[], double b[])
/**C*F****************************************************************
**                                                                  **
** Assigns D-dimensional vector b to vector a.                      **
** You might encounter problems with the macro ASSIGND on some      **
** machines. If yes, better use this function although it's slower. **
**                                                                  **
***C*F*E*************************************************************/
{
   int j;
   for (j=0; j<D; j++)
   {
      a[j] = b[j];
   }
}

// 產生 0 ~ 1 間的亂數
double rnd_uni(long *idum)
/**C*F****************************************************************
**                                                                  **
** SRC-FUNCTION   :rnd_uni()                                        **
** LONG_NAME      :random_uniform                                   **
** AUTHOR         :(see below)                                      **
**                                                                  **
** DESCRIPTION    :rnd_uni() generates an equally distributed ran-  **
**                 dom number in the interval [0,1]. For further    **
**                 reference see Press, W.H. et alii, Numerical     **
**                 Recipes in C, Cambridge University Press, 1992.  **
**                                                                  **
** FUNCTIONS      :none                                             **
**                                                                  **
** GLOBALS        :none                                             **
**                                                                  **
** PARAMETERS     :*idum    serves as a seed value                  **
**                                                                  **
** PRECONDITIONS  :*idum must be negative on the first call.        **
**                                                                  **
** POSTCONDITIONS :*idum will be changed                            **
**                                                                  **
***C*F*E*************************************************************/
{
  long j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  double temp;

  if (*idum <= 0)
  {
    if (-(*idum) < 1) *idum=1;
    else *idum = -(*idum);
    idum2=(*idum);
    for (j=NTAB+7;j>=0;j--)
    {
      k=(*idum)/IQ1;
      *idum=IA1*(*idum-k*IQ1)-k*IR1;
      if (*idum < 0) *idum += IM1;
      if (j < NTAB) iv[j] = *idum;
    }
    iy=iv[0];
  }
  k=(*idum)/IQ1;
  *idum=IA1*(*idum-k*IQ1)-k*IR1;
  if (*idum < 0) *idum += IM1;
  k=idum2/IQ2;
  idum2=IA2*(idum2-k*IQ2)-k*IR2;
  if (idum2 < 0) idum2 += IM2;
  j=iy/NDIV;
  iy=iv[j]-idum2;
  iv[j] = *idum;
  if (iy < 1) iy += IMM1;
  if ((temp=AM*iy) > RNMX) return RNMX;
  else return temp;

}/*------End of rnd_uni()--------------------------*/

// 將上下限轉為全域變數
double inibound_h;      /* upper parameter bound              */
double inibound_l;      /* lower parameter bound              */
// 與機構合成相關的全域變數
// 宣告一個座標結構
struct Coord {
    double x;
    double y;
  // 這裡保留 double z;
};

main(int argc, char *argv[])
/**C*F****************************************************************
**                                                                  **
** SRC-FUNCTION   :main()                                           **
** LONG_NAME      :main program                                     **
** AUTHOR         :Rainer Storn, Kenneth Price                      **
**                                                                  **
** DESCRIPTION    :driver program for differential evolution.       **
**                                                                  **
** FUNCTIONS      :rnd_uni(), evaluate(), printf(), fprintf(),      **
**                 fopen(), fclose(), fscanf().                     **
**                                                                  **
** GLOBALS        :rnd_uni_init    input variable for rnd_uni()     **
**                                                                  **
** PARAMETERS     :argc            #arguments = 3                   **
**                 argv            pointer to argument strings      **
**                                                                  **
** PRECONDITIONS  :main must be called with three parameters        **
**                 e.g. like de1 <input-file> <output-file>, if     **
**                 the executable file is called de1.               **
**                 The input file must contain valid inputs accor-  **
**                 ding to the fscanf() section of main().          **
**                                                                  **
** POSTCONDITIONS :main() produces consecutive console outputs and  **
**                 writes the final results in an output file if    **
**                 the program terminates without an error.         **
**                                                                  **
***C*F*E*************************************************************/

{
   char  chr;             /* y/n choice variable                */
   char  *strat[] =       /* strategy-indicator                 */
   {
            "",
            "DE/best/1/exp",
            "DE/rand/1/exp",
            "DE/rand-to-best/1/exp",
            "DE/best/2/exp",
            "DE/rand/2/exp",
            "DE/best/1/bin",
            "DE/rand/1/bin",
            "DE/rand-to-best/1/bin",
            "DE/best/2/bin",
            "DE/rand/2/bin"
   };

   int   i, j, L, n;      /* counting variables                 */
   int   r1, r2, r3, r4;  /* placeholders for random indexes    */
   int   r5;              /* placeholders for random indexes    */
   int   D;               /* Dimension of parameter vector      */
   int   NP;              /* number of population members       */
   int   imin;            /* index to member with lowest energy */
   int   refresh;         /* refresh rate of screen output      */
   int   strategy;        /* choice parameter for screen output */
   int   gen, genmax, seed;   

   long  nfeval;          /* number of function evaluations     */

   double trial_cost;      /* buffer variable                    */
   // 將上下限轉為全域變數, 可能要根據各變數加以設定
   //double inibound_h;      /* upper parameter bound              */
   //double inibound_l;      /* lower parameter bound              */
   double tmp[MAXDIM], best[MAXDIM], bestit[MAXDIM]; /* members  */
   double cost[MAXPOP];    /* obj. funct. values                 */
   double cvar;            /* computes the cost variance         */
   double cmean;           /* mean cost                          */
   double F,CR;            /* control variables of DE            */
   double cmin;            /* help variables                     */

   FILE  *fpin_ptr;
   FILE  *fpout_ptr;

// 計算執行過程所需時間起點, 需要導入 time.h
  clock_t start = clock();

/*------Initializations----------------------------*/

 //if (argc != 3)                                 /* number of arguments */
 //{
    //printf("\nUsage : de <input-file> <output-file>\n");
    //exit(1);
 //}

// 將結果寫入 out.dat
 fpout_ptr = fopen("out.dat","w");          /* open output file for reading,    */
                                          /* to see whether it already exists */
 /*
 if ( fpout_ptr != NULL )
 {
    printf("\nOutput file %s does already exist, \ntype y if you ",argv[2]);
    printf("want to overwrite it, \nanything else if you want to exit.\n");
    chr = (char)getchar();
    if ((chr != 'y') && (chr != 'Y'))
    {
      exit(1);
    }
    fclose(fpout_ptr);
 }
*/

/*-----Read input data------------------------------------------------*/

 //fpin_ptr   = fopen(argv[1],"r");
/*
 if (fpin_ptr == NULL)
 {
    printf("\nCannot open input file\n");
    exit(1);
 }*/

 //fscanf(fpin_ptr,"%d",&strategy);       /*---choice of strategy-----------------*/
 //fscanf(fpin_ptr,"%d",&genmax);         /*---maximum number of generations------*/
 //fscanf(fpin_ptr,"%d",&refresh);        /*---output refresh cycle---------------*/
 //fscanf(fpin_ptr,"%d",&D);              /*---number of parameters---------------*/
 //fscanf(fpin_ptr,"%d",&NP);             /*---population size.-------------------*/
 //fscanf(fpin_ptr,"%lf",&inibound_h);    /*---upper parameter bound for init-----*/
 //fscanf(fpin_ptr,"%lf",&inibound_l);    /*---lower parameter bound for init-----*/
 //fscanf(fpin_ptr,"%lf",&F);             /*---weight factor----------------------*/
 //fscanf(fpin_ptr,"%lf",&CR);            /*---crossing over factor---------------*/
 //fscanf(fpin_ptr,"%d",&seed);           /*---random seed------------------------*/
// 目前已經採用 strategy 3 可以得到最佳結果
  strategy = 3;
  genmax = 2000;
  refresh = 100;
  // 配合機構尺寸合成, 每一個體有 9 個機構尺寸值與 5 個通過點角度值
  D = 19;
  NP = 200;
  inibound_h = 50.;
  inibound_l = 0.;
/*得到最佳解
  F = 0.85;
CR 必須介於 0 to 1. 之間
  CR = 1.;
*/
  F = 0.85;
  CR = 1.;
  seed = 3;

 //fclose(fpin_ptr);

/*-----Checking input variables for proper range----------------------------*/

  if (D > MAXDIM)
  {
     printf("\nError! D=%d > MAXDIM=%d\n",D,MAXDIM);
     exit(1);
  }
  if (D <= 0)
  {
     printf("\nError! D=%d, should be > 0\n",D);
     exit(1);
  }
  if (NP > MAXPOP)
  {
     printf("\nError! NP=%d > MAXPOP=%d\n",NP,MAXPOP);
     exit(1);
  }
  if (NP <= 0)
  {
     printf("\nError! NP=%d, should be > 0\n",NP);
     exit(1);
  }
  if ((CR < 0) || (CR > 1.0))
  {
     printf("\nError! CR=%f, should be ex [0,1]\n",CR);
     exit(1);
  }
  if (seed <= 0)
  {
     printf("\nError! seed=%d, should be > 0\n",seed);
     exit(1);
  }
  if (refresh <= 0)
  {
     printf("\nError! refresh=%d, should be > 0\n",refresh);
     exit(1);
  }
  if (genmax <= 0)
  {
     printf("\nError! genmax=%d, should be > 0\n",genmax);
     exit(1);
  }
  if ((strategy < 0) || (strategy > 10))
  {
     printf("\nError! strategy=%d, should be ex {1,2,3,4,5,6,7,8,9,10}\n",strategy);
     exit(1);
  }
  if (inibound_h < inibound_l)
  {
     printf("\nError! inibound_h=%f < inibound_l=%f\n",inibound_h, inibound_l);
     exit(1);
  }

/*-----Open output file-----------------------------------------------*/

   //fpout_ptr   = fopen(argv[2],"w");  /* open output file for writing */

   //if (fpout_ptr == NULL)
   //{
      //printf("\nCannot open output file\n");
      //exit(1);
   //}

/*-----Initialize random number generator-----------------------------*/

 rnd_uni_init = -(long)seed;  /* initialization of rnd_uni() */
 nfeval       =  0;  /* reset number of function evaluations */

/*------Initialization------------------------------------------------*/
/*------Right now this part is kept fairly simple and just generates--*/
/*------random numbers in the range [-initfac, +initfac]. You might---*/
/*------want to extend the init part such that you can initialize-----*/
/*------each parameter separately.------------------------------------*/

   for (i=0; i<NP; i++)
   {
      for (j=0; j<D; j++) /* spread initial population members */
      {
        c[i][j] = inibound_l + rnd_uni(&rnd_uni_init)*(inibound_h - inibound_l);
      }
      cost[i] = evaluate(D,c[i],&nfeval); /* obj. funct. value */
   }
   cmin = cost[0];
   imin = 0;
   for (i=1; i<NP; i++)
   {
     if(MAXIMAPROBLEM == 1)
     {
       // 改為最大化
        if (cost[i]>cmin)
        {
          cmin = cost[i];
          imin = i;
        }
      }
      else
      {
        // 最小化問題
        if (cost[i]<cmin)
        {
          cmin = cost[i];
          imin = i;
        }
      }
   }

   assignd(D,best,c[imin]);            /* save best member ever          */
   assignd(D,bestit,c[imin]);          /* save best member of generation */

   pold = &c; /* old population (generation G)   */
   pnew = &d; /* new population (generation G+1) */

/*=======================================================================*/
/*=========Iteration loop================================================*/
/*=======================================================================*/

   gen = 0;                          /* generation counter reset */
   while ((gen < genmax) /*&& (kbhit() == 0)*/) /* remove comments if conio.h */
   {                                            /* is accepted by compiler    */
      gen++;
      imin = 0;

      for (i=0; i<NP; i++)         /* Start of loop through ensemble  */
      {
     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 2 !!!     */
       r1 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while(r1==i);            

     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 3 !!!     */
       r2 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while((r2==i) || (r2==r1));

     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 4 !!!     */
       r3 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while((r3==i) || (r3==r1) || (r3==r2));

     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 5 !!!     */
       r4 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while((r4==i) || (r4==r1) || (r4==r2) || (r4==r3));

     do                        /* Pick a random population member */
     {                         /* Endless loop for NP < 6 !!!     */
       r5 = (int)(rnd_uni(&rnd_uni_init)*NP);
     }while((r5==i) || (r5==r1) || (r5==r2) || (r5==r3) || (r5==r4));

/*=======Choice of strategy===============================================================*/
/*=======We have tried to come up with a sensible naming-convention: DE/x/y/z=============*/
/*=======DE :  stands for Differential Evolution==========================================*/
/*=======x  :  a string which denotes the vector to be perturbed==========================*/
/*=======y  :  number of difference vectors taken for perturbation of x===================*/
/*=======z  :  crossover method (exp = exponential, bin = binomial)=======================*/
/*                                                                                        */
/*=======There are some simple rules which are worth following:===========================*/
/*=======1)  F is usually between 0.5 and 1 (in rare cases > 1)===========================*/
/*=======2)  CR is between 0 and 1 with 0., 0.3, 0.7 and 1. being worth to be tried first=*/
/*=======3)  To start off NP = 10*D is a reasonable choice. Increase NP if misconvergence=*/
/*           happens.                                                                     */
/*=======4)  If you increase NP, F usually has to be decreased============================*/
/*=======5)  When the DE/best... schemes fail DE/rand... usually works and vice versa=====*/

/*=======EXPONENTIAL CROSSOVER============================================================*/

/*-------DE/best/1/exp--------------------------------------------------------------------*/
/*-------Our oldest strategy but still not bad. However, we have found several------------*/
/*-------optimization problems where misconvergence occurs.-------------------------------*/
     if (strategy == 1) /* strategy DE0 (not in our paper) */
     {
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D);
       L = 0;
       do
       {                       
         tmp[n] = bestit[n] + F*((*pold)[r2][n]-(*pold)[r3][n]);
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }
/*-------DE/rand/1/exp-------------------------------------------------------------------*/
/*-------This is one of my favourite strategies. It works especially well when the-------*/
/*-------"bestit[]"-schemes experience misconvergence. Try e.g. F=0.7 and CR=0.5---------*/
/*-------as a first guess.---------------------------------------------------------------*/
     else if (strategy == 2) /* strategy DE1 in the techreport */
     {
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D);
       L = 0;
       do
       {                       
         tmp[n] = (*pold)[r1][n] + F*((*pold)[r2][n]-(*pold)[r3][n]);
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }
/*-------DE/rand-to-best/1/exp-----------------------------------------------------------*/
/*-------This strategy seems to be one of the best strategies. Try F=0.85 and CR=1.------*/
/*-------If you get misconvergence try to increase NP. If this doesn't help you----------*/
/*-------should play around with all three control variables.----------------------------*/
     else if (strategy == 3) /* similiar to DE2 but generally better */
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
       L = 0;
       do
       {                       
         tmp[n] = tmp[n] + F*(bestit[n] - tmp[n]) + F*((*pold)[r1][n]-(*pold)[r2][n]);
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }
/*-------DE/best/2/exp is another powerful strategy worth trying--------------------------*/
     else if (strategy == 4)
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
       L = 0;
       do
       {                           
         tmp[n] = bestit[n] + 
              ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F;
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }
/*-------DE/rand/2/exp seems to be a robust optimizer for many functions-------------------*/
     else if (strategy == 5)
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
       L = 0;
       do
       {                           
         tmp[n] = (*pold)[r5][n] + 
              ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F;
         n = (n+1)%D;
         L++;
       }while((rnd_uni(&rnd_uni_init) < CR) && (L < D));
     }

/*=======Essentially same strategies but BINOMIAL CROSSOVER===============================*/

/*-------DE/best/1/bin--------------------------------------------------------------------*/
     else if (strategy == 6) 
     {
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = bestit[n] + F*((*pold)[r2][n]-(*pold)[r3][n]);
         }
         n = (n+1)%D;
           }
     }
/*-------DE/rand/1/bin-------------------------------------------------------------------*/
     else if (strategy == 7) 
     {
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = (*pold)[r1][n] + F*((*pold)[r2][n]-(*pold)[r3][n]);
         }
         n = (n+1)%D;
           }
     }
/*-------DE/rand-to-best/1/bin-----------------------------------------------------------*/
     else if (strategy == 8) 
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = tmp[n] + F*(bestit[n] - tmp[n]) + F*((*pold)[r1][n]-(*pold)[r2][n]);
         }
         n = (n+1)%D;
           }
     }
/*-------DE/best/2/bin--------------------------------------------------------------------*/
     else if (strategy == 9)
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = bestit[n] + 
              ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F;
         }
         n = (n+1)%D;
           }
     }
/*-------DE/rand/2/bin--------------------------------------------------------------------*/
     else
     { 
       assignd(D,tmp,(*pold)[i]);
       n = (int)(rnd_uni(&rnd_uni_init)*D); 
           for (L=0; L<D; L++) /* perform D binomial trials */
           {
         if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */
         {                       
           tmp[n] = (*pold)[r5][n] + 
              ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F;
         }
         n = (n+1)%D;
           }
     }

/*=======Trial mutation now in tmp[]. Test how good this choice really was.==================*/

     trial_cost = evaluate(D,tmp,&nfeval);  /* Evaluate new vector in tmp[] */
   if(MAXIMAPROBLEM == 1)
   {
    // 改為最大化
       if (trial_cost >= cost[i])   /* improved objective function value ? */
       {                                  
          cost[i]=trial_cost;         
          assignd(D,(*pnew)[i],tmp);
          if (trial_cost>cmin)          /* Was this a new minimum? */
          {                               /* if so...*/
             cmin=trial_cost;           /* reset cmin to new low...*/
             imin=i;
             assignd(D,best,tmp);           
          }                           
       }                            
       else
       {
          assignd(D,(*pnew)[i],(*pold)[i]); /* replace target with old value */
       }
    }
    else
    {
          // 最小化問題
       if (trial_cost <= cost[i])   /* improved objective function value ? */
       {                                  
          cost[i]=trial_cost;         
          assignd(D,(*pnew)[i],tmp);
          if (trial_cost<cmin)          /* Was this a new minimum? */
          {                               /* if so...*/
             cmin=trial_cost;           /* reset cmin to new low...*/
             imin=i;
             assignd(D,best,tmp);           
          }                           
       }                            
       else
       {
          assignd(D,(*pnew)[i],(*pold)[i]); /* replace target with old value */
       }
    }

      }   /* End mutation loop through pop. */

      assignd(D,bestit,best);  /* Save best population member of current iteration */

      /* swap population arrays. New generation becomes old one */

      pswap = pold;
      pold  = pnew;
      pnew  = pswap;

/*----Compute the energy variance (just for monitoring purposes)-----------*/

      cmean = 0.;          /* compute the mean value first */
      for (j=0; j<NP; j++)
      {
         cmean += cost[j];
      }
      cmean = cmean/NP;

      cvar = 0.;           /* now the variance              */
      for (j=0; j<NP; j++)
      {
         cvar += (cost[j] - cmean)*(cost[j] - cmean);
      }
      cvar = cvar/(NP-1);

/*----Output part----------------------------------------------------------*/

      if (gen%refresh==1)   /* display after every refresh generations */
      { /* ABORT works only if conio.h is accepted by your compiler */
    printf("\n\n                         PRESS ANY KEY TO ABORT"); 
    printf("\n\n\n Best-so-far cost funct. value=%-15.10g\n",cmin);

    for (j=0;j<D;j++)
    {
      printf("\n best[%d]=%-15.10g",j,best[j]);
    }
    printf("\n\n Generation=%d  NFEs=%ld   Strategy: %s    ",gen,nfeval,strat[strategy]);
    printf("\n NP=%d    F=%-4.2g    CR=%-4.2g   cost-variance=%-10.5g\n",
               NP,F,CR,cvar);
      }

      fprintf(fpout_ptr,"%ld   %-15.10g\n",nfeval,cmin);
   }
/*=======================================================================*/
/*=========End of iteration loop=========================================*/
/*=======================================================================*/

/*-------Final output in file-------------------------------------------*/

   fprintf(fpout_ptr,"\n\n\n Best-so-far obj. funct. value = %-15.10g\n",cmin);

   for (j=0;j<D;j++)
   {
     fprintf(fpout_ptr,"\n best[%d]=%-15.10g",j,best[j]);
   }
   fprintf(fpout_ptr,"\n\n Generation=%d  NFEs=%ld   Strategy: %s    ",gen,nfeval,strat[strategy]);
   fprintf(fpout_ptr,"\n NP=%d    F=%-4.2g    CR=%-4.2g    cost-variance=%-10.5g\n",
           NP,F,CR,cvar); 

  fclose(fpout_ptr);

  /* Code you want timed here */
  printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
   return(0);
}

/*-----------End of main()------------------------------------------*/

// 適應函式 fittness function (cost function)
double evaluate(int D, double tmp[], long *nfeval)
{
  // 先處理通過 5 個點的四連桿問題
  // x0, y0 為左方固定點座標, 必須在 0 ~ 100 間 - 設為 tmp[0], tmp[1]
  // x1, y1 為右方固定點座標, 必須在 0 ~ 100 間 - 設為 tmp[2], tmp[3]
  // L1 為第一桿件的長度, 必須 > 0, 且小於 100 - 設為 tmp[4]
  // L2 為第二桿件的長度, 必須 > 0, 且小於 100 - 設為 tmp[5]
  // L3 為第三桿件的長度, 必須 > 0, 且小於 100 - 設為 tmp[6]
  // L5, L6 與 L2 共同圍出可動桿 L2 對應的三角形, 關注的三角形頂點即 L5 與 L6 的交點, 而 angle3 則為 L6 之對應角(為固定值)
  // L5, L6 必須 > 0, 且小於 100 - 設為 tmp[7], tmp[8]
  // 以下的角度輸入值, 會隨著目標點數的增加而增加, 其索引值由 9 + 通過點數 - 1 決定, 5 點, 則索引至 13, 若通過 25 點, 則索引值為 9 + 24 = 33
  // input_angles[] 為五的輸入的雙浮點數角度值,代表個體的角度向量值 - 分別設為 tmp[9], tmp[10], tmp[11], tmp[12], tmp[13]
  // output_points[] 則為與 input_angles[] 對應的五個三角形的頂點值, 為座標結構值, 分別有 x 與 y 分量值
  // 當利用個體的向量值, 代入 mechanism 後所得到得 output_points[] 再與 target_points[] 進行 cost function 的誤差值最小化
  /* void mechanism(double x0, double y0, double x1, double y1, double L1,
  double L2, double L3, double L5, double L6, double input_angles[NUM_OF_POINTS], struct Coord output_points[NUM_OF_POINTS])*/
  struct Coord target_points[NUM_OF_POINTS], output_points[NUM_OF_POINTS];
  double input_angles[NUM_OF_POINTS], result;
  int i;

  (*nfeval)++;

  target_points[0].x = 1.0;
  target_points[0].y = 1.0;

  target_points[1].x = 2.0;
  target_points[1].y = 2.0;

  target_points[2].x = 3.0;
  target_points[2].y = 3.0;

  target_points[3].x = 4.0;
  target_points[3].y = 4.0;

  target_points[4].x = 5.0;
  target_points[4].y = 5.0;

  target_points[5].x = 6.0;
  target_points[5].y = 6.0;

  target_points[6].x = 7.0;
  target_points[6].y = 7.0;

  target_points[7].x = 8.0;
  target_points[7].y = 8.0;

  target_points[8].x = 9.0;
  target_points[8].y = 9.0;

  target_points[9].x = 10.0;
  target_points[9].y = 10.0;

  // 輸入角度值與 tmp[] 的設定
  for(i = 0; i < NUM_OF_POINTS; i++)
  {
    input_angles[i] = tmp[i + 9];
  }
  // 呼叫 mechanism() 以便計算 output_points[]
  mechanism(tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], tmp[6], tmp[7], tmp[8], input_angles, output_points);

  // for debug
  /*
  if(*nfeval%3000 == 0)
  {
    for(i = 0; i < NUM_OF_POINTS; i++)
    {
      printf("%-15.10g : %-15.10g\n", output_points[i].x, output_points[i].y);
    }
    printf("#####################################\n");
  }
  */
  // 利用 output_points[] 與 target_points 計算誤差值, 該誤差值就是 cost
  result = error_function(output_points, target_points);
  // 這裡要分別針對各變數的約束條件與範圍值來設定傳回 PENALITY 或 誤差值 result

  // x0 與 x1 點位於 -500 與 500 中間
    for(i = 0; i < 4; i++)
  {
    if(tmp[i] < -50 || tmp[i] > 50){
      return PENALITY;
    }
  }

  // 三個連桿值, 一定要為正
    for(i = 4; i < 7; i++)
  {
    if(tmp[i] < 0 || tmp[i] > 50){
      return PENALITY;
    }
  }

    // L5 L6 可以為 0 或負值
    for(i = 7; i < 9; i++)
  {
    if(tmp[i] < -50 || tmp[i] > 50){
      return PENALITY;
    }
  }

  // 角度值一定要大於 0

  for(i = 1; i <= NUM_OF_POINTS; i++)
  {
    if((tmp[D-i] < 0)){
      return PENALITY;
    }
  }

  return result;

  /*
   double result=0, surface = 80.0, z, volume, penality;
   (*nfeval)++;
   z = (surface-tmp[0]*tmp[1])/(2.0*(tmp[0]+tmp[1]));
   volume = tmp[0]*tmp[1]*z;

  if(volume <= 0){
    return PENALITY;
  }

  if((tmp[0] <= inibound_l)|| (tmp[0] >inibound_h)){
    return PENALITY;
  }

  if((tmp[1] <= inibound_l) || (tmp[1] >inibound_h)){
    return PENALITY;
  }
  // volume must >0 and max volume
  // 目前為最小化問題
   return 1+1/(volume*volume);
   */
}

struct Coord triangletip_coord( double x0, double y0, double R0, double R1, double x1, double y1, double localt)
{
    struct Coord tip_coord;

    if (localt>=0 && localt <PI)
    {
        // 目前蓋掉的式子為利用手動代換出來的版本
        //x_value = ((x1-x0)*(-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2))/(2*sqrt(pow((y1-y0),2)+pow((x1-x0),2)))-(y1-y0)*(-mech_loop)*sqrt(fabs(pow(R0,2)-pow((-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2)),2)/(4*(pow((y1-y0),2)+pow((x1-x0),2))))))/sqrt(pow((y1-y0),2)+pow((x1-x0),2))+x0;
        // 以下的式子,先利用文字編輯器,將原先 stringout() 出來的  sqrt 替換成  sqrtt, 以防止被  maxima 中的 subst("^"=pow,expr) 所替換, subst 之後,再使用文字編輯器換回來,就可以得到正確的 C 對應碼.
        tip_coord.x = pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*(mech_loop*(y1-y0)*sqrt(fabs(pow(R0,2)-pow(pow(y1-y0,2)+
    pow(x1-x0,2),-1)*pow(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2),2)/4))+(x1-x0)*pow(sqrt(pow(y1-y0,2)+
    pow(x1-x0,2)),-1)*(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2))/2)+x0;
    }
    else
    {
        // 目前蓋掉的式子為利用手動代換出來的版本
        //x_value = ((x1-x0)*(-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2))/(2*sqrt(pow((y1-y0),2)+pow((x1-x0),2)))-(y1-y0)*(mech_loop)*sqrt(fabs(pow(R0,2)-pow((-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2)),2)/(4*(pow((y1-y0),2)+pow((x1-x0),2))))))/sqrt(pow((y1-y0),2)+pow((x1-x0),2))+x0;
        tip_coord.x = pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*(-mech_loop*(y1-y0)*sqrt(fabs(pow(R0,2)-pow(pow(y1-y0,2)+
    pow(x1-x0,2),-1)*pow(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2),2)/4))+(x1-x0)*pow(sqrt(pow(y1-y0,2)+
    pow(x1-x0,2)),-1)*(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+
    pow(x1-x0,2))/2)+x0;
    }

// 請注意,與 Maxma 公式中的差異為,在 sqrt()中加入 fabs(),避免因為sqrt()中的負值而造成 NaN (Not a Number 問題.
    if (localt>=0 && localt <PI)
    {
        tip_coord.y = /*((x1-x0)*(-mech_loop)*sqrt(
                fabs(pow(R0,2)-pow((-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2)),2)
                /(4*(pow((y1-y0),2)+pow((x1-x0),2)))
                ))
                +(y1-y0)*(-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2))/(2*sqrt(pow((y1-y0),2)+pow((x1-x0),2))))/sqrt(pow((y1-y0),2)+pow((x1-x0),2))
                +y0;*/
                // 利用 sqrtt 居中進行代換所得到的式子
                pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*((y1-y0)*pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*(-pow(R1,2)+
    pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2))/2-mech_loop*(x1-x0)*sqrt(fabs(pow(R0,2)-pow(pow(y1-y0,2)+
    pow(x1-x0,2),-1)*pow(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2),2)/4)))+y0;

    }
    else
    {
        tip_coord.y = /*((x1-x0)*(mech_loop)*sqrt(
                fabs(pow(R0,2)-pow((-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2)),2)
                /(4*(pow((y1-y0),2)+pow((x1-x0),2)))
                ))
                +(y1-y0)*(-pow(R1,2)+pow(R0,2)+pow((y1-y0),2)+pow((x1-x0),2))/(2*sqrt(pow((y1-y0),2)+pow((x1-x0),2))))/sqrt(pow((y1-y0),2)+pow((x1-x0),2))
                +y0;*/
                pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*((y1-y0)*pow(sqrt(pow(y1-y0,2)+pow(x1-x0,2)),-1)*(-pow(R1,2)+
    pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2))/2+mech_loop*(x1-x0)*sqrt(fabs(pow(R0,2)-pow(pow(y1-y0,2)+
    pow(x1-x0,2),-1)*pow(-pow(R1,2)+pow(R0,2)+pow(y1-y0,2)+pow(x1-x0,2),2)/4)))+y0;
    }

  return tip_coord;
}

double distance(double x0, double y0, double x1, double y1)
{
    double distance_value;
    distance_value = sqrt(pow((x1-x0),2) + pow((y1-y0),2));
    return distance_value;
}

double rr(double L1, double dd, double theta)
{
    double rr_value;
    rr_value = sqrt(L1*L1+dd*dd-2*L1*dd*cos(theta));
    return rr_value;
}

// 輸入每一個體的變數向量, 然後求各三角形頂點的座標陣列[NUM_OF_POINTS]
void mechanism(double x0, double y0, double x1, double y1, double L1,
  double L2, double L3, double L5, double L6, double input_angles[NUM_OF_POINTS], struct Coord output_points[NUM_OF_POINTS])
{
  // 此函式要輸入控制變數, 然後計算機構尺寸合成的關注點座標位置
  // 以下為可能的處理變數宣告
  // 這裡希望能夠定義一個 struct 來處理座標點
  double rr_length, dd_length, angle;
  struct Coord link1_tip, link2_tip, triangle_tip;
    double angle2, angle3;
  int i;

  // 開始進行三角形頂點座標的計算
  // 以下變數由每一個體向量提供
  /*
    x0 = 0.0;
    y0 = 0.0;
    x1 = 10.0;
    y1 = 0.0;
    L1 = 5.0;
    L2 = 20;
    L3 = 10;
    L5 = 10;
    L6 = 10;
  */
  dd_length = distance(x0, y0, x1, y1);
  /* 設法表示 triangle 所對應的 local 角度,表示為已知變數與 t 的函式 */
  angle3 = acos((pow(L2,2)+pow(L5,2)-pow(L6,2))/(2*L2*L5));

  for(i = 0; i < NUM_OF_POINTS; i++)
  {
    // 先建立第一點座標, 即 i=0 者
    // i=0;
    // angle = i*degree;
    /*
    // 利用角度增量進行運算, 相對於 input_angles[0] 作為基準
    if(i > 0)
    {
      input_angles[i] = input_angles[i] + input_angles[i-1];
    }
    */
    angle = input_angles[i]*degree;
    rr_length = rr(L1, dd_length, angle);
    // 第一次三角形疊代
    link1_tip = triangletip_coord(x0, y0, L1, rr_length, x1, y1, angle);
    // 第二次三角形疊代
    /* 設法表示 link2 所對應的 local 角度,表示為已知變數與 t 的函式 */
    angle2 = acos((pow(L2,2)+pow(rr_length,2)-pow(L3,2))/(2*L2*rr_length));
    link2_tip = triangletip_coord(link1_tip.x, link1_tip.y, L2, L3, x1, y1, angle2);
    // 第三次三角形疊代
    //triangle_tip = triangletip_coord(link1_tip.x, link1_tip.y, L5, L6, link2_tip.x, link2_tip.y, angle3);
    // output_points[i] = triangletip_coord(link1_tip.x, link1_tip.y, L5, L6, link2_tip.x, link2_tip.y, angle3);
    // 這裡要嘗試利用 finaltip_coord() 求 tip3 座標, 而 L5 與 L6 可 0 可負
    output_points[i] = finaltip_coord(link1_tip, link2_tip, L5, L6);
  }
}

double error_function(struct Coord output_points[NUM_OF_POINTS], struct Coord target_points[NUM_OF_POINTS])
{
  double error = 0.0;
  int i;
  for(i = 0; i < NUM_OF_POINTS; i++)
  {
    error += fabs(distance(output_points[i].x, output_points[i].y, target_points[i].x, target_points[i].y));
  }
  return error;
}

struct Coord finaltip_coord(struct Coord tip1_coord, struct Coord tip2_coord, double r1, double r2)
{
  struct Coord tip3_coord;
  double theta3, theta4, length3, length4;
  length3 = sqrt(pow(tip2_coord.x - tip1_coord.x,2) + pow(tip2_coord.y - tip1_coord.y,2));
  length4 = sqrt(pow(r1,2) + pow(r2,2));  
  theta3 = acos((tip2_coord.x - tip1_coord.x) / length3);
  theta4 = acos(r1 / length4);
  tip3_coord.x = tip1_coord.x + length4 * cos(theta3 + theta4);
  tip3_coord.y = tip1_coord.y + length4 * sin(theta3 + theta4);

  return tip3_coord;
}