/*
  このプログラムをtype1.cpp"という名前で保存して、
  gcc -g v-ebata.cpp -o v-ebata
  でコンパイルして下さい。

===============================================================

    ダイエット仮想人格シミュレータ「バーチャル江端」
   (ファジィ推論エンジン + 超シンプル体重変動シミュレータ拡張版)

                                                    2016/02/07
                                                    江端智一

    このプログラムの利用によって発生する不利益対して、
    江端の完全な免責を条件として、このプログラムの
    無償かつ無制限の使用を許諾します。
*/

#include <stdio.h>
#include <stdlib.h>

#define LAST14   14  // 14日分の数字"14"を使う時

// 共通関数
double max(double a){
  return a;
}

double max(double a, double b){
  if (a > b) 
    return a;
  else 
    return b;
};

double max(double a, double b, double c ){
  return max(max(a,b), max(b,c)); 
};

double max(double a, double b, double c, double d ){
  return max(max(a,b,c), d); 
};

double min(double a){
  return a;
}

double min(double a, double b){
  if (b > a) 
    return a;
  else 
    return b;
};

double min(double a, double b, double c ){
  return min(min(a,b), min(b,c)); 
};

double min(double a, double b, double c, double d ){
  return min(min(a,b,c), d); 
};


// ファジィ表現
typedef enum scale {LESSLESS, LESS, ZERO, MORE, MOREMORE} SCALE;

// 前件部メンバーシップ関数(山3つ)クラス
class condition_MF3
{
private:
  double center;
  double width;
  SCALE express;
  
public:
  condition_MF3(double _center, double _witdth, SCALE _express){
    center = _center;
    width = _witdth;
    express = _express;
    
    // 使用できないファジィ表現を使った場合は止める        
    if ((express == LESSLESS) || (express == MOREMORE)){
      printf("wrong expression used \n");
      exit(0);
    }
    
  };
  double func(double _x);
};

double condition_MF3::func(double _x)
{
  // x,yは、メンバーシップ関数上の座標を示す
  double x = _x;
  double y = 0.0; // yの値は、必ず0以上1以下になる
  
  if (express == LESS){
    if (x <= center - width){
      y = 1.0;
    }
    else if (x <= center){
      y = - 1.0 / width * (x - center);
    }
    else{
      y = 0.0;
    }
  }
  else if (express == ZERO){
    if (x <= center - width){
      y = 0.0;
    }
    else if (x <= center){
      y = 1.0 / width * (x - center) + 1.0;
    }
    else if (x <= center + width){
      y = -1.0 / width * (x - center) + 1.0;
    }
    else{
      y = 0.0;
    }
  }
  else if (express == MORE){
    if (x <= center){
      y = 0.0;
    }
    else if (x <= center + width){
      y = 1.0 / width * (x - center);
    }
    else{
      y = 1.0;
    }
  }
  else {
    printf("wrong expression\n");
    exit(1);
  }
  
  return y;
};

// 前件部メンバーシップ関数(山5つ)クラス
class condition_MF5
{
private:
  double center;
  double width;
  SCALE express;
  
public:
  condition_MF5(double _center, double _witdth, SCALE _express){
    center = _center;
    width = _witdth;
    express = _express;
  };
  double func(double _x);
};


double condition_MF5::func(double _x)
{
  // x,yは、メンバーシップ関数上の座標を示す
  double x = _x;
  double y = 0.0; // yの値は、必ず0以上1以下になる
  
  if (express == LESSLESS){
    if (x <= center - 2.0 * width){
      y = 1.0;
    }
    else if (x <= center - width){
      y = - 1.0 / width * (x - (center - 2.0 * width)) + 1.0;
    }
    else{
      y = 0.0;
    }
  }
  else if (express == LESS){
    if (x <= center - 2.0 * width){
      y = 0.0;
    }
    else if (x <= center - width){
      y = 1.0 / width * (x - (center - width)) + 1.0;
    }
    else if (x <= center){
      y = -1.0 / width * (x - (center - width)) + 1.0; 
    }
    else{
      y = 0.0;
    }
  }
  else if (express == ZERO){
    if (x <= center - width){
      y = 0.0;
    }
    else if (x <= center){
      y = 1.0 / width * (x - center) + 1.0;
    }
    else if (x <= center + width){
      y = -1.0 / width * (x - center) + 1.0;
    }
    else{
      y = 0.0;
    }
  }
  else if (express == MORE){
    if (x <= center){
      y = 0.0;
    }
    else if (x <= center + width){
      y = 1.0 / width * (x - (center + width)) + 1.0;
    }
    else if (x <= center + 2.0 * width){
      y = -1.0 / width * (x - (center + width)) + 1.0; 
    }
    else{
      y = 0.0;
    }
  }
  else if (express == MOREMORE){
    if (x <= center + width){
      y = 0.0;
    }
    else if (x <= center + 2.0 * width){
      y = 1.0 / width * (x - (center + 2.0 * width)) + 1.0;
    }
    else{
      y = 1.0;
    }
  }
  
  return y;
};

// 後件部メンバーシップ関数(山3つ)クラス  
class action_MF3
{
private:
  double center;
  double width;
  SCALE express;
  
  double x;
  double y;
  
public:
  action_MF3(double _center, double _witdth, SCALE _express){
    
    y = 0.0; // yの値は、必ず0以上1以下になる
    
    center = _center;
    width = _witdth;
    express = _express;
    
    if (express == LESS){
      x = center - width;
    }
    else if (express == ZERO){
      x = center;
    }
    else if (express == MORE){
      x = center + width;
    }
    else{
      printf("wrong scale expression\n");
      exit(0);
    }
  };

  // Y座標の値を最大値で更新する  
  void func_Max(double b){
    y = max(b, y);
  };
  
  // Y座標の値をリセット(y=0)する
  void func_Reset(){
    y = 0.0;
  };
  
  // X座標を返す
  double func_X(void){
    return x;
  };
  
  // (最大値で更新された、最後の)Y座標を返す
  double func_Y(){
    return y;
  };

};

// 後件部メンバーシップ関数(山5つ)クラス  
class action_MF5
{
private:
  double center;
  double width;
  SCALE express;
  
  double x;
  double y;
  
public:
  action_MF5(double _center, double _witdth, SCALE _express){
    y = 0.0; // yの値は、必ず0以上1以下になる
    
    center = _center;
    width = _witdth;
    express = _express;
    
    if (express == LESSLESS){
      x = center - 2.0 * width;
    }
    else if (express == LESS){
      x = center - width;
    }
    else if (express == ZERO){
      x = center;
    }
    else if (express == MORE){
      x = center + width;
    }
    else if (express == MOREMORE){
      x = center + 2.0 * width;
    }
    else{
      printf("wrong scale expression\n");
      exit(-1); // 強制終了
    }
  };
  
  // Y座標の値を最大値で更新する  
  void func_Max(double b){
    y = max(b, y);
  };
  
  // Y座標の値をリセット(y=0)する
  void func_Reset(){
    y = 0.0;
  };
  
  // X座標を返す
  double func_X(void){
    return x;
  };
  
  // (最大値で更新された、最後の)Y座標を返す
  double func_Y(){
    return y;
  };
  
};

int main()
{
  // ダイエット経過期間
  condition_MF3 Period_Short(14, 14, LESS); // 0日
  condition_MF3 Period_Middle(14, 14, ZERO); // 14日
  condition_MF3 Period_Long(14, 14, MORE); // 28日
  
  // いわゆる「停滞期」(14日間あたりの平均減重変化量)
  condition_MF3 WeighChange_Small(1.0, 1.0, LESS); // 0kg
  condition_MF3 WeighChange_Middle(1.0, 1.0, ZERO); // 1kg
  condition_MF3 WeighChange_High(1.0, 1.0, MORE); // 2kg
  
  // 苦痛度
  condition_MF3 Suffer_Zerol(400.0, 400.0, LESS); // 0
  condition_MF3 Suffer_Middle(400.0, 400.0, ZERO); // 400
  condition_MF3 Suffer_High(400.0, 400.0, MORE); // 800
  
  // BMI
  condition_MF5 BMI_lowlow(22.0, 3.0, LESSLESS); // 痩せすぎ(15.0)
  condition_MF5 BMI_low(22.0, 3.0, LESS);// 痩せぎみ(19.0)
  condition_MF5 BMI_Norcmal(22.0, 3.0, ZERO);// 普通(22.0)
  condition_MF5 BMI_High(22.0, 3.0, MORE); // 太りぎみ(25.0)
  condition_MF5 BMI_HighHigh(22.0, 3.0, MOREMORE);// 太りすぎ(28.0)
  
  // 一日の摂取カロリー
  action_MF5 Eat_LittleLittle(2500, 1000, LESSLESS); // 500 Kcal
  action_MF5 Eat_Little(2500, 1000, LESS);// 1500 Kcal
  action_MF5 Eat_Normal(2500, 1000, ZERO);// 2500 Kcal
  action_MF5 Eat_Lot(2500, 1000, MORE);// 3500 Kcal
  action_MF5 Eat_LotLot(2500, 1000, MOREMORE);// 4500 Kcal
  
  // 超シンプルシミュレータ投入
  // 体重変動シミュレーションプログラム
  //基礎代謝量と求め方 ?ハリス・ベネディクト方程式(日本人版)　計算式
  //http://www.kintore.info/kisotaisya_mass/index.cgi
  //【計算式】
  // 男性 66.5+(体重kg×13.8)+(身長cm×5.0)-(年齢×6.8)
  // 女性 665.1+(体重kg×9.6)+(身長cm×1.9)-(年齢×7.0)                         
  // 江端智一 2014日2月4日現在のデータを使って計算
  double weight = 78.0;  // 開始時の体重 78.0kg
  double lengthw = 172.0;  // 身長 172cm
  double age = 49.0 + 94/365;  // 開始時の年齢 49歳と94日 (2014年2月4日現在)
  
  // 最初の基礎代謝カロリー 男性  66.5+(体重kg×13.8)+(身長cm×5.0)-(年齢×6.8)
  double consumption_calorie =  
    66.5 + weight * 13.8 + 172.0 * 5.0 - age * 6.8;
  // このカロリーに運動消費カロリーを加えて、
  // トータルの一日の消費カロリーを算出する
  
  // (現在の体重 66.5kgにおいて)一日運動消費カロリー 708kcalと推定し、
  // このカロリーは、体重に比例するものと仮定する。
  consumption_calorie += weight/66.5 * 708.0;
  //printf("consumption_calorie = %f\n",  consumption_calorie);
  
  // 摂取カロリーの初期値は、暫定的に消費カロリーと同値とする
  // (第一回目の苦痛値の算出の為)
  double source_calorie = consumption_calorie;
  
  double last_weight[LAST14]; // 過去14日分の体重
  
  // 最初は全部同じ体重とする 14日経過までは同じ体重とする
  for (int i = 0; i < LAST14; i++ ){
    last_weight[i] = weight;
  }
  
  printf("経過日,体重,消費カロリー,摂取カロリー\n");
  
  // ここからループ開始
  for (int day=0; day < 30; day++){ // ここで日数を入れてね(ここでは30日)
    // 7kcal = 1g = 0.001kg とする    
    
    // 体重データを一日分ずらす
    for (int i = 0;  i < LAST14-1; i++ ){
      last_weight[i+1] = last_weight[i];
    }
    last_weight[0] = weight;
    
    // 線形補完の傾きだけを求める(14日分)
    // http://d.hatena.ne.jp/rainlib/20090112/1231735459
    double k1=0, k2=0, k3=0, k4=0;
    
    for (int i = 0; i < LAST14; i++ ){
      k1 += (double)i * (double)i;         // Σ(x^2)
      k2 += (double)i * last_weight[i];    // Σ(xy)
      k3 += (double)i;                     // Σ(x) 
      k4 += last_weight[i];                // Σ(y)
    }
    
    double dW = ((double)LAST14 * k2 - k3 * k4) / ((double)LAST14 * k4 - k3*k3); // 14
    
    // BMIの計算式 体重(Kg)÷身長(m)÷身長(m) 江端の身長172cm
    double bmi = weight / 1.72 / 1.72;
    
    // 苦痛値の算出
    double suffer = consumption_calorie - source_calorie;
    //printf("suffer=%f\n", suffer);
    
    double r; // min-maxのmin値算出結果用の変数
    
    
    // =============== ファジイルール入れ替え ===============

    // =============== バーチャル江端 タイプ1(ここから) ===============
   
    // [ルール1] ダイエット開始直後は、「無条件にがんばれる」
    r = min(Period_Short.func(day));
    Eat_Little.func_Max(r);
    //printf("r1: %f\n",r);
    
    // [ルール2] しかし、苦しいと感じたら、食べ過ぎてしまう。
    r = min(Suffer_Middle.func(suffer));
    Eat_Lot.func_Max(r);
    //printf("r2: %f\n",r);    
    
    r = min(Suffer_High.func(suffer));
    Eat_Lot.func_Max(r);
    //printf("r2: %f\n",r);    
    
    // [ルール3] 後半になると、ダイエットのことを忘れて食べてしまう
    r = min(Period_Long.func(day));
    Eat_Lot.func_Max(r);
    //printf("r3: %f\n",r);    
    
    // =============== バーチャル江端 タイプ1(ここまで) ===============
    
    /*
      ルールに振れない場合は、min-max重心法の分母がゼロになり、
      ゼロ割が発生する場合がある為、それを回避する
    */
    double denominator =  // 分母
      (Eat_LittleLittle.func_Y() + 
       Eat_Little.func_Y() +
       Eat_Normal.func_Y() +
       Eat_Lot.func_Y() +
       Eat_LotLot.func_Y()) ;
    
    double numerator =  // 分子
      (Eat_LittleLittle.func_X() * Eat_LittleLittle.func_Y() + 
       Eat_Little.func_X() * Eat_Little.func_Y() +
       Eat_Normal.func_X() * Eat_Normal.func_Y() +
       Eat_Lot.func_X() * Eat_Lot.func_Y() +
       Eat_LotLot.func_X() * Eat_LotLot.func_Y());
    
    // 推論結果 (分母がゼロの場合は、推論結果は前回と同値とする)
    if ( denominator != 0.0){
      source_calorie =  numerator / denominator ;
    }      
    
    // 後件部メンバーシップ関数のリセット(ループさせる時は必ずリセットする)
    
    Eat_LittleLittle.func_Reset();
    Eat_Little.func_Reset();
    Eat_Normal.func_Reset();
    Eat_Lot.func_Reset();
    Eat_LotLot.func_Reset();
    
    //printf("\n source_calorie = %f\n",source_calorie);
    
    double diff_weight = (source_calorie - consumption_calorie)/7.0 / 1000.0;
    weight += diff_weight;
    age += 1.0/365.0;
    consumption_calorie =  66.5 + weight * 13.8 + 172.0 * 5.0 - age * 6.8;
    consumption_calorie += weight/66.5 * 708.0;
    
    // printf("day:%d\tweight = %f\tconsumption_calorie =%f\n", day,weight, consumption_calorie);
    
    printf("%d\t%5.2f\t%5.2f\t%5.2f\n", day,weight, consumption_calorie,source_calorie);
  }
}

/*
  必要に応じてルールを差し替えて下さい。

    // =============== バーチャル江端 タイプ1(ここから) ===============
   
    // [ルール1] ダイエット開始直後は、「無条件にがんばれる」
    r = min(Period_Short.func(day));
    Eat_Little.func_Max(r);
    //printf("r1: %f\n",r);
    
    // [ルール2] しかし、苦しいと感じたら、食べ過ぎてしまう。
    r = min(Suffer_Middle.func(suffer));
    Eat_Lot.func_Max(r);
    //printf("r2: %f\n",r);    
    
    r = min(Suffer_High.func(suffer));
    Eat_Lot.func_Max(r);
    //printf("r2: %f\n",r);    
    
    // [ルール3] 後半になると、ダイエットのことを忘れて食べてしまう
    r = min(Period_Long.func(day));
    Eat_Lot.func_Max(r);
    //printf("r3: %f\n",r);    
    
    // =============== バーチャル江端 タイプ1(ここまで) ===============

    // =============== バーチャル江端 タイプ2(ここから) ===============

    // [ルール1] 1月程度は「がんばれる」
    r = min(Period_Short.func(day));
    Eat_Little.func_Max(r);
    //printf("r1: %f\n",r);

    r = min(Period_Middle.func(day));
    Eat_Little.func_Max(r);
    //printf("r1: %f\n",r);

    
    // [ルール2] しかし、苦しいと感じた日には普通に食べてしまうこともある。
    r = min(Suffer_Middle.func(suffer));
    Eat_Normal.func_Max(r);
    //printf("r2: %f\n",r);    

    // [ルール3] 停滞が続くと食べてしまう。
    r = min(WeighChange_Middle.func(dW));
    Eat_Normal.func_Max(r);

    // [ルール4] 全く停滞してしまうと、食べてしまう。
    r = min(WeighChange_Small.func(dW));
    Eat_Normal.func_Max(r);
    
    // =============== バーチャル江端 タイプ2(ここまで) ===============

    // =============== バーチャル江端 タイプ3(ここから) ===============

    // [ルール1] いつだってがんばっている
    r = min(Period_Short.func(day));
    Eat_Little.func_Max(r);
    //printf("r1: %f\n",r);

    r = min(Period_Middle.func(day));
    Eat_Little.func_Max(r);
    //printf("r1: %f\n",r);

    r = min(Period_Middle.func(day));
    Eat_Little.func_Max(r);
    //printf("r1: %f\n",r);

    // [ルール2] やせてくると、コントロールできなくなる
    r = min(BMI_low.func(bmi));
    Eat_LittleLittle.func_Max(r);

    r = min(BMI_lowlow.func(bmi));
    Eat_LittleLittle.func_Max(r);

    // =============== バーチャル江端 タイプ3(ここまで) ===============
  
*/