プログラマ専用SNS ミクプラ

コード:

#define _CRT_RAND_S

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <math.h>
#include "VariationalBayes.h"
#include <float.h>
#include <windows.h>
#include "gamma.h"
#include "utility.h"
#include <shlobj.h>

#pragma warning(disable: 4996)


// 0～1の乱数を返す
double RandF()
{
	unsigned int rnd;

	return (double)rnd/UINT_MAX;
}

// 整数の乱数を返す
unsigned int RandD()
{
	unsigned int rnd;
	rnd = rand();
	return rnd;
}

CVariationalBayes::CVariationalBayes()
{
	//構造体初期化
	memset( m_models , 0 , MAXMODEL*sizeof(MODEL) );

	for(int i=0;i<DIM;i++)
	{
		m_data[i] = NULL;
	}
}


CVariationalBayes::~CVariationalBayes()
{
	FreeMemory();
}



// プログラム全体を通じて最後に一度実行される関数
void CVariationalBayes::FreeMemory()
{
	for(int i=0 ; i<MAXMODEL ; i++) FreeModel( &(m_models[i]) );
}

// モデル構造体を確保
int CVariationalBayes::AllocModel(MODEL *model , int modelnum )
{
	model->m = new PARAMETER[modelnum+1];

	for(int j=0 ; j<=modelnum ;j++)
	{
		model->m[j].z = new double[m_numData];
	}

	model->m_freeEnergies = 0;
	model->ModelNum = modelnum+1;

	return 0;
}

int CVariationalBayes::FreeModel( MODEL *model )
{
	if( model->m )
	{
		for(int j=0 ; j<model->ModelNum ;j++)
		{
			delete [] model->m[j].z;
			model->m[j].z = NULL;
		}

		delete [] model->m;
		model->m = NULL;
	}

	model->m_freeEnergies = 0;
	model->ModelNum = 0;

	return 0;
}

// 学習関数
void CVariationalBayes::Learn( std::string SaveDir )
{
	//データ保存用ディレクトリ
	if( SaveDir.size() )
	{
		if( *(SaveDir.end()-1) != '/' || *(SaveDir.end()-1) != '\\' ) SaveDir += "/";
	}
	else
	{
		printf("ディレクトリを指定してください。");
		return;
	}
	CreateDirectory( SaveDir.c_str() , NULL );


	//メモリ確保
	for(int i=0 ; i<MAXMODEL ; i++)
	{
		if( m_models[i].m ) FreeModel( &(m_models[i]) );
		AllocModel( &(m_models[i]) , i );
		m_freeEnergies[i] = -DBL_MAX;
	}


	//model数を変えて最適事後分布を計算
	for(int M=1;M<MAXMODEL;M++)
	{
		printf("モデル数 %02d 計算開始" , M+1 );

		PARAMETER *params = m_models[M].m;

		for(int it=0 ; it<ITERATE ; it++ )
		{
			double oldF = -DBL_MAX;
			double Fm = 0;
			
			//Step1-1:初期パラメータ決定
			InitializeParameter(M);

			//Step1-2:変数初期値決定
			for(int m=0;m<=M;m++){
				SetInitializingParameter(&params[m],M,m);
			}

			for(int loop=0 ; ; loop++ )
			{
				//Step2-1:潜在変数の変分事後分布の更新
				for(int n=0;n<m_numData;n++)
				{
					double gamma[MAXMODEL];

					// γの計算
					for(int m=0;m<=M;m++)
					{
						gamma[m]=CalcGamma(M,m,n);
					}

					// \bar{z}の計算
					for(int a=0;a<=M;a++)
					{
						double sum=0;
						for(int k=0;k<=M;k++)
						{
							sum+=exp(gamma[k]-gamma[a]);
						}
						params[a].z[n]=1.0/sum;
					}
				}

				//Step2-2:パラメータの変分事後分布の更新
				for(int m=0;m<=M;m++){
					PARAMETER &p = params[m];
					double invBeta[DIM][DIM] = {0};

					//初期化
					p.Nbar=0;
					for(int l=0;l<DIM;l++){
						p.xbar[l]=0;
						for(int ll=0;ll<DIM;ll++)
							p.C[l][ll]=0;
					}

					// \bar{x}_i, \bar{N}_i の計算
					for(int n=0;n<m_numData;n++){
						p.Nbar+=p.z[n];
						for(int l=0;l<DIM;l++)
							p.xbar[l]+=p.z[n]*m_data[l][n];
					}
					for(int b=0;b<DIM;b++)
						p.xbar[b]/=p.Nbar;

					// \bar{C}_i の計算
					for(int c=0;c<m_numData;c++){
						for(int l=0;l<DIM;l++)
							for(int ll=0;ll<DIM;ll++){
								p.C[l][ll] += p.z[c]*(m_data[l][c]-p.xbar[l])*(m_data[ll][c]-p.xbar[ll]);
							}
					}

					// φ, η, fの更新
					p.phi=m_phi0+p.Nbar;
					p.eta=m_eta0+p.Nbar;
					p.f=p.eta+1-DIM;

					//混合比
					p.alpha=p.Nbar/m_numData;

					// \bar{μ}, β, Σの更新
					InverseMatrix(p.beta,invBeta,DIM);	//+++
					for(int d=0;d<DIM;d++){
						p.mubar[d]=(p.Nbar*p.xbar[d]+m_xi0*m_nu0[m][d])/(p.Nbar+m_xi0);
						p.mu[d]=p.mubar[d];	//平均
						for(int ll=0;ll<DIM;ll++){
							p.beta[d][ll]=m_beta0[d][ll]+p.C[d][ll]+(p.Nbar*m_xi0)/(p.Nbar+m_xi0)*(p.xbar[d]-m_nu0[m][d])*(p.xbar[ll]-m_nu0[m][ll]);
							p.sigma[d][ll]=p.beta[d][ll]/((p.Nbar+m_xi0)*p.f);
							p.ivariance[d][ll]=p.eta*invBeta[d][ll];		//分散の逆行列(精度行列)(ウィシャート分布の期待値)
						}
					}
					InverseMatrix(p.ivariance,p.variance,DIM);			//分散(variance)行列
					for(int e=0;e<DIM;e++)
						for(int ll=0;ll<DIM;ll++)
							p.deviation[e][ll]=sqrt(p.variance[e][ll]);	//標準偏差行列
				}

				Fm = CalcFreeEnergy( m_models[M] , M );

				// 収束判定
				if( Fm - oldF < 0.00000001 )
				{
					break;
				}

				oldF = Fm;
			}

			// 自由エネルギーが最大の物を選択
			if( m_freeEnergies[M] < Fm )
			{
				m_freeEnergies[M] = Fm;
				SaveModelData( SaveDir , m_models[M],m_freeEnergies[M],M);	// 計算結果保存
			}
		}

		printf("計算完了...F = %lf\n" , m_freeEnergies[M] );
	}

	// 最適なモデルの選択
	double maxF = -DBL_MAX;
	for(int A=1 ; A<MAXMODEL ; A++ )
	{
		if( m_freeEnergies[A] > maxF )
		{
			maxF = m_freeEnergies[A];
			m_bestModelIdx = A;
		}
	}

}


// データをセット
int CVariationalBayes::SetData(double **data , int numData )
{
	//メモリ確保
	m_numData = numData;


	for(int i=0;i<DIM;i++)
	{
		if( m_data[i] ) free( m_data[i] );

		m_data[i] = new double[m_numData];
	}


	//データをコピー
	for(int cnt=0 ; cnt<m_numData ; cnt++)
	{
		for(int i=0;i<DIM;i++)
		{
			m_data[i][cnt] = data[cnt][i];
		}
	}

	return 0;
}


//初期パラメータの値決め
int CVariationalBayes::InitializeParameter(int model )
{
	m_phi0 = m_numData;
	m_xi0 = 1;
	m_eta0 = DIM + 2;

	double mu[DIM] = {0};	// 各次元の平均
	double var[DIM] = {0};	// 各次元の分散

	// 各次元の平均
	for(int n=0 ; n<m_numData ; n++ )
		for(int d=0 ; d<DIM ; d++ )	mu[d] += m_data[d][n];
	for(int d=0 ; d<DIM ; d++ )	mu[d] /= m_numData;

	// 各次元の分散
	for(int f=0 ; f<m_numData ; f++ )
		for(int d=0 ; d<DIM ; d++ )	var[d] += (m_data[d][f] - mu[d])*(m_data[d][f] - mu[d]);
	for(int g=0 ; g<DIM ; g++ )	var[g] /= m_numData;


	// ランダムでデータを選択
	for(int m=0;m<=model;m++)
	{
		int rnd = RandD() % m_numData;
		for(int d=0; d<DIM ; d++)
			m_nu0[m][d] = m_data[d][rnd];
	}


	// 分散
	for(int d1=0; d1<DIM ; d1++)
		for(int d0=0; d0<DIM ; d0++)
		{
			if( d0==d1 ) m_beta0[d1][d0] = var[d0];
			else m_beta0[d1][d0] = 0;
		}


	return 0;
}



//初期パラメータを代入する
void CVariationalBayes::SetInitializingParameter( PARAMETER* prm, int maxmodel, int model )
{
	prm->Nbar=(double)m_numData/(maxmodel+1);
	prm->phi=m_phi0 + prm->Nbar;
	prm->eta=m_eta0 + prm->Nbar;
	prm->f=m_eta0+prm->Nbar+1-DIM;

	for(int l=0;l<DIM;l++)
	{
		prm->mubar[l]=m_nu0[model][l];
		for(int ll=0;ll<DIM;ll++)
		{
			prm->beta[l][ll]=m_beta0[l][ll];
			prm->sigma[l][ll]=prm->beta[l][ll]/((prm->Nbar+m_xi0)*prm->f);
		}
	}
	return;
}

//γ[m][m_numData]を計算する
double CVariationalBayes::CalcGamma( int maxmodel, int model, int dataNo)
{
	double p1,p2,p3,p4,p5_1[DIM][DIM],p5_2[DIM][DIM],sum_p5[DIM][DIM],p5;
	double sunNbar,sumDigamma,iBeta[DIM][DIM],track;

	sunNbar = 0;
	sumDigamma = 0;

	for(int i=0;i<=maxmodel;i++)
		sunNbar+=m_models[maxmodel].m[i].Nbar;

	p1=digamma(m_phi0+m_models[maxmodel].m[model].Nbar);
	p2=-1*digamma((maxmodel+1)*m_phi0+sunNbar);

	for(int d=1;d<=DIM;d++)
		sumDigamma+=digamma((m_models[maxmodel].m[model].eta-d)/2);

	p3=0.5*sumDigamma;
	p4=(-0.5)*log(Determinant(m_models[maxmodel].m[model].beta,DIM));
	InverseMatrix(m_models[maxmodel].m[model].beta,iBeta,DIM);

	for(int h=0;h<DIM;h++)
	{
		for(int j=0;j<DIM;j++)
		{
			sum_p5[h][j]=0;
			p5_1[h][j]=m_models[maxmodel].m[model].eta*iBeta[h][j];
			p5_2[h][j]=(m_models[maxmodel].m[model].f/(m_models[maxmodel].m[model].f-2.0)*m_models[maxmodel].m[model].sigma[h][j])+(m_data[h][dataNo]-m_models[maxmodel].m[model].mubar[h])*(m_data[j][dataNo]-m_models[maxmodel].m[model].mubar[j]);
		}
	}

	for(int l=0;l<DIM;l++)
		for(int ll=0;ll<DIM;ll++)
			for(int k=0;k<DIM;k++)
				sum_p5[l][ll]+=p5_1[l][k]*p5_2[k][ll];

	track = TraceMatrix(sum_p5,DIM);
	p5 = (-0.5)*track;

	return p1+p2+p3+p4+p5;
}


//|A|:行列式
double CVariationalBayes::Determinant(double A[][DIM], int dim)
{
	double t,u,det;
	double S[DIM][DIM];

	for(int i=0;i<dim;i++)
		for(int j=0;j<dim;j++)
			S[i][j]=A[i][j];

	det = 1;
	for(int k=0;k<dim;k++){
		t = S[k][k];
		det *=t;
		for(int i=0;i<dim;i++) S[k][i] /=t;
		S[k][k] = 1/t;
		for(int j=0;j<dim;j++)
			if(j!=k){
				u=S[j][k];
				for(int i=0;i<dim;i++)
					if(i!=k) S[j][i] -= S[k][i] * u;
					else	 S[j][i] = -u/t;
			}
	}
	return det;
}

//逆行列:入力行列→A　出力用逆行列→S
void CVariationalBayes::InverseMatrix(double A[][DIM],double S[][DIM], int dim)
{
	double t,u,det;

	for(int i=0;i<dim;i++)
		for(int j=0;j<dim;j++)
			S[i][j]=A[i][j];

	det = 1;
	for(int k=0;k<dim;k++){
		t = S[k][k];
		det *=t;
		for(int i=0;i<dim;i++) S[k][i] /=t;
		S[k][k] = 1/t;
		for(int j=0;j<dim;j++)
			if(j!=k){
				u=S[j][k];
				for(int i=0;i<dim;i++)
					if(i!=k) S[j][i] -= S[k][i] * u;
					else	 S[j][i] = -u/t;
			}
	}
}

//行列のトレース::対角成分の和を計算
double CVariationalBayes::TraceMatrix(double S[][DIM], int dim)
{
	double sum=0;
	for(int i=0;i<dim;i++)
		sum+=S[i][i];
	return sum;
}

//自由エネルギーを計算する
double CVariationalBayes::CalcFreeEnergy( MODEL model, int m)
{
	double term[10] = {0};
	PARAMETER *p = model.m;
	double xi[MAXMODEL] = {0};
	double detB[MAXMODEL] = {0};
	double F = 0;

	// \xi_iの計算
	for(int i=0 ; i<=m ; i++ ) xi[i] = m_xi0 + p[i].Nbar;

	// |B_0|, |B_i|の計算
	double detB0 = Determinant( m_beta0 , DIM );
	for(int k=0 ; k<=m ; k++ ) detB[k] = Determinant( p[k].beta , DIM );

	// 第一項目
	for(int o=0 ; o<=m ; o++ )  term[0] += gamma( log( p[o].phi ) );

	// 第二項目
	term[1] = gamma( log( (m+1) * m_phi0) );

	// 第三項目
	term[2] = -gamma( log( (m+1)*m_phi0 + m_numData ) );

	// 第四項目
	term[3] = -gamma( log(m_phi0) ) * (m+1);

	// 第五項目
	term[4] = -( DIM * m_numData * log(m_phi0)/log(2.0) )/2.0;

	//term[4] = -( DIM * m_numData * log(M_PI)/log(2.0) )/2.0;

	// 第六項目
	for(int r=0 ; r<=m ; r++ ) term[5] += log( m_xi0 / xi[r] );
	term[5] *= DIM/2.0;

	// 第七項目
	term[6] = m_eta0 * log( detB0 ) / 2.0;
	for(int D=1 ; D<=DIM ; D++ ) term[6] -= gamma( log( (m_eta0+1-D) /2.0) );
	term[6] *= (m+1);

	// 第八項目
	for(int s=0 ; s<=m ; s++ )
	{
		term[7] += p[s].eta * log(detB[s]) /2.0;
		for(int D=1 ; D<=DIM ; D++ )
		{
			term[7] -= gamma( log( (p[s].eta+1-D)/2.0 ) );
		}
	}
	term[7] *= -1;

	// 総和
	for(int t=0 ; t<8 ; t++ )
	{
		F += term[t];
	}

	return F;
}

// 結果をファイルに保存する
void CVariationalBayes::SaveModelData( std::string saveDir , MODEL mod, double F , int modNum )
{
	double zmax,zresult;
	char strNo[32];
	char strBuff[256];
	sprintf( strNo , "%02d" , modNum+1 );
	FILE *fp = fopen( (saveDir + "result" + strNo + ".txt").c_str() , "w" );
	std::string filename = saveDir + "detail" + strNo + ".txt";


	// モデルのパラメータの保存
	SaveVariable( "■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■" , "設定値" ,  filename.c_str() , "w");
	SaveVariable( modNum+1 , "モデル数" , filename.c_str() );
	SaveVariable( m_numData , "データ数" , filename.c_str() );
	SaveVariable( F , "自由エネルギー" , filename.c_str() );

	SaveVariable( "■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■" , "\n\n初期値" , filename.c_str() );
	SaveVariable( m_phi0 , "φ" , filename.c_str() );
	SaveVariable( m_eta0 , "η" , filename.c_str() );
	SaveVariable( m_xi0 , "ξ" , filename.c_str() );
	SaveVariableMatrix( m_nu0 , modNum+1 , DIM , "ν" , filename.c_str() );

	for(int m=0;m<=modNum;m++)
	{
		sprintf( strBuff , "\n\nmodel = %d" , m );
		SaveVariable( "■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■" , strBuff , filename.c_str() );

		sprintf( strBuff , "Nbar[%d]" , m );
		SaveVariable( mod.m[m].Nbar , strBuff , filename.c_str() );

		sprintf( strBuff , "phi[%d]" , m );
		SaveVariable( mod.m[m].phi , strBuff , filename.c_str() );

		sprintf( strBuff , "eta[%d]" , m );
		SaveVariable( mod.m[m].eta, strBuff , filename.c_str() );

		sprintf( strBuff , "f[%d]" , m );
		SaveVariable( mod.m[m].f , strBuff , filename.c_str() );

		sprintf( strBuff , "mubar[%d]" , m );
		SaveVariableArray( mod.m[m].mubar , DIM , strBuff , filename.c_str() );

		sprintf( strBuff , "xbar[%d]" , m );
		SaveVariableArray( mod.m[m].xbar , DIM , strBuff , filename.c_str() );

		sprintf( strBuff , "Beta[%d]" , m );
		SaveVariableMatrix( mod.m[m].beta , DIM , DIM , strBuff , filename.c_str() );

		sprintf( strBuff , "Sigma[%d]" , m );
		SaveVariableMatrix( mod.m[m].sigma , DIM , DIM , strBuff , filename.c_str() );

		sprintf( strBuff , "Cbar[%d]" , m );
		SaveVariableMatrix( mod.m[m].C , DIM , DIM , strBuff , filename.c_str() );

		sprintf( strBuff , "\nα[%d]" , m );
		SaveVariable( mod.m[m].alpha , strBuff , filename.c_str() );

		sprintf( strBuff , "平均[%d]" , m );
		SaveVariableArray( mod.m[m].mu , DIM , strBuff , filename.c_str() );

		sprintf( strBuff , "共分散[%d]" , m );
		SaveVariableMatrix( mod.m[m].variance , DIM , DIM , strBuff , filename.c_str() );

		sprintf( strBuff , "逆共分散[%d]" , m );
		SaveVariableMatrix( mod.m[m].ivariance , DIM , DIM , strBuff , filename.c_str() );
	}

	// クラスタリング結果を保存
	for(int n=0;n<m_numData;n++)
	{
		zresult=zmax=0;
		for(int m=0;m<=modNum;m++)
		{
			fprintf( fp ,"%.5lf	",mod.m[m].z[n]);
			if(zmax<mod.m[m].z[n])
			{
				zmax=mod.m[m].z[n];	
				zresult=m+1;	
			}
		}
		fprintf( fp ,"%.1lf\n",zresult);
	}
	fclose( fp );

}

コード:

#pragma once
/*******************************************************************

utility.h		developed by naka_t	2011.01.31

	ファイル入出力系関数

　＊VariationalBayes.hに必要な関数を移植			naka_t	2011.01.31

Copyright (C) 2011  naka_t <naka_t@apple.ee.uec.ac.jp>
*******************************************************************/
#include <vector>
#include <windows.h>
#include <iostream>

#pragma warning(disable: 4996)
#pragma warning(disable: 4018)

#ifndef Message
#define Message( mess ) ::MessageBox( 0 , mess , "Message" , 0 )
#endif


// ファイル入出力関数
static int fprintfVar( int val , FILE *fp ){ return fprintf( fp , "%d" , val ); }
static int fprintfVar( double val , FILE *fp ){ return fprintf( fp , "%lf" , val ); }
static int fprintfVar( std::string &val , FILE *fp ){ return fprintf( fp , "%s" , val.c_str() ); }
static int fprintfVar( const char *&val , FILE *fp ){ return fprintf( fp , "%s" , val ); }
static int fscanfVar( int *var , FILE *fp ){ return fscanf( fp, "%d" , var ); }
static int fscanfVar( double *var , FILE *fp ){ return fscanf( fp, "%lf" , var ); }
static int fscanfVar( float *var , FILE *fp ){ return fscanf( fp, "%f" , var ); }
static int fscanfVar( bool *var , FILE *fp ){ return fscanf( fp, "%d" , var ); }
static int fscanfVar( std::string *var , FILE *fp ){
	char tmpStr[1024] = {0};
	int ret = fscanf( fp , "%s" , tmpStr );
	*var = tmpStr;
	return ret;
}

// std::string版のfgets
// 戻り値：読み込んだ文字数
//		   ファイルの終端の際は、string::nposを返します
static size_t fgets( std::string &buff , FILE *fp ){
	char c;
	buff.clear();
	while ( ( c = fgetc( fp ) ) != EOF && c != '\n' )
		buff.push_back( c );

	if( c==EOF && buff.size()==0 ) return std::string::npos;

	return buff.size();
}

//名前つきで変数を保存
template <typename T> bool SaveVariable( T value , const char *valname , const char *filename , const char *mode="a" )
{
	FILE *fp;
	std::vector<std::string> vecLines;
	std::string line;
	bool bFind = false;

	if( *mode=='a')	// 追記モードならば現在の内容を一時保存
	{
		fp = fopen( filename , "r" );
		if( fp )	// ファイルがあれば読み込む
		{
			while( fgets( line , fp ) != std::string::npos ) { vecLines.push_back(line); }
			fclose(fp);
		}
	}


	// 書き込み開始
	fp = fopen( filename , "w" );

	if( !fp )
	{
		char mess[1024];
		sprintf( mess , "SaveVariable : %sがオープンできません" , filename );
		Message( mess );
		return false;
	}

	for(int i=0 ; i<vecLines.size() ; i++ )
	{
		char tmpName[256];
		sscanf( vecLines[i].c_str() , "%s" , tmpName );
		if( strcmp(tmpName,valname)==0 )
		{
			// 変数が見つかったら、上書き
			fprintf( fp , "%s		" , valname );
			fprintfVar( value , fp );
			fprintf( fp , "\n" );
			bFind = true;
		}
		else
		{
			// 見つからなかったら、そのまま保存
			fprintf( fp , "%s\n" , vecLines[i].c_str() );
		}
	}

	if( !bFind )
	{
		fprintf( fp , "%s		" , valname );
		fprintfVar( value , fp );
		fprintf( fp , "\n" );
	}

	fclose(fp);
	return true;
}


//名前つきで変数を保存
template <typename T> bool SaveVariableArray( T value , int num , const char *valname , const char *filename , const char *mode="a" )
{
	FILE *fp;
	std::vector<std::string> vecLines;
	std::string line;
	bool bFind = false;

	if( *mode=='a')	// 追記モードならば現在の内容を一時保存
	{
		fp = fopen( filename , "r" );
		if( fp )	// ファイルがあれば読み込む
		{
			while( fgets( line , fp ) != std::string::npos ) { vecLines.push_back(line); }
			fclose(fp);
		}
	}


	// 書き込み開始
	fp = fopen( filename , "w" );

	if( !fp )
	{
		char mess[1024];
		sprintf( mess , "SaveVariable : %sがオープンできません" , filename );
		Message( mess );
		return false;
	}

	for(int i=0 ; i<vecLines.size() ; i++ )
	{
		char tmpName[256];
		sscanf( vecLines[i].c_str() , "%s" , tmpName );
		if( strcmp(tmpName,valname)==0 )
		{
			// 変数が見つかったら、上書き
			fprintf( fp , "%s	" , valname );
			for(int j=0 ; j<num ; j++ )
			{
				fprintfVar( value[j] , fp );
				fprintf( fp , "	" );
			}
			fprintf( fp , "\n" );
			bFind = true;
		}
		else
		{
			// 見つからなかったら、そのまま保存
			fprintf( fp , "%s\n" , vecLines[i].c_str() );
		}
	}

	if( !bFind )
	{
		fprintf( fp , "%s	" , valname );
		for(int j=0 ; j<num ; j++ )
		{
			fprintfVar( value[j] , fp );
			fprintf( fp , "	" );
		}
		fprintf( fp , "\n" );
	}

	fclose(fp);
	return true;
}

//名前つきで変数を保存
template <typename T> bool SaveVariableMatrix( T value , int ysize , int xsize , const char *valname , const char *filename , const char *mode="a" )
{
	FILE *fp;
	std::vector<std::string> vecLines;
	std::string line;
	bool bFind = false;

	if( *mode=='a')	// 追記モードならば現在の内容を一時保存
	{
		fp = fopen( filename , "r" );
		if( fp )	// ファイルがあれば読み込む
		{
			while( fgets( line , fp ) != std::string::npos ) { vecLines.push_back(line); }
			fclose(fp);
		}
	}


	// 書き込み開始
	fp = fopen( filename , "w" );

	if( !fp )
	{
		char mess[1024];
		sprintf( mess , "SaveVariable : %sがオープンできません" , filename );
		Message( mess );
		return false;
	}

	for(int i=0 ; i<vecLines.size() ; i++ )
	{
		char tmpName[256];
		sscanf( vecLines[i].c_str() , "%s" , tmpName );
		if( strcmp(tmpName,valname)==0 )
		{
			fprintf( fp , "%s			%d	%d\n" , valname , ysize , xsize );

			// 変数が見つかったら、上書き
			for(int y=0 ; y<ysize ; y++ )
			{
				fprintf( fp , "*	" );
				for(int x=0 ; x<xsize ; x++ )
				{
					fprintfVar( value[y][x] , fp );
					fprintf( fp , "	" );
				}
				fprintf( fp , "\n" );
			}

			// "*"が行列の続きを表すので読み飛ばす
			for( i = i+1 ; i<vecLines.size() ; i++ )
			{
				sscanf( vecLines[i].c_str() , "%s" , tmpName );
				if( strcmp(tmpName,"*")!=0 ) 
				{
					i--;
					break;
				}
			}
			bFind = true;
		}
		else
		{
			// 見つからなかったら、そのまま保存
			fprintf( fp , "%s\n" , vecLines[i].c_str() );
		}
	}

	if( !bFind )
	{
		fprintf( fp , "%s			%d	%d\n" , valname , ysize , xsize );
		for(int y=0 ; y<ysize ; y++ )
		{
			fprintf( fp , "*	" );
			for(int x=0 ; x<xsize ; x++ )
			{
				fprintfVar( value[y][x] , fp );
				fprintf( fp , "	" );
			}
			fprintf( fp , "\n" );
		}

	}

	fclose(fp);
	return true;
}


//行数をカウント
static int CountLine( FILE *fp )
{
	int num = 0;
	int c;

	while(1)
	{
		c = fgetc(fp);

		if( c == '\0' || c == EOF ) break;
		else if( c == '\n' ) num++;
	}

	//printf( "lines %d\n" , num );

	rewind(fp);

	return num;
}

//行内のタブの数を数える
static int CountCols( FILE *fp )
{
	int num = 0;
	int c , old;

	while(1)
	{
		c = fgetc(fp);

		if( c == '\n' || c == EOF ) break;
		else if( c == '\t' ) num++;

		old = c;
	}

	//最後の列にタブが入っていない可能性もあるので
	if( old != '\t' && c == '\n' ) num++;

	rewind(fp);

	return num;
}

// メモリ確保
template<typename T> T** AllocMatrix( int ysiz , int xsiz)
{
	//配列の実態は1次元で確保していることに注意！

	T **mptr = NULL;
	T *data = NULL;	//配列の実態

	try
	{
		mptr = new T*[ysiz];
		data = new T[ysiz*xsiz];

		for (int i=0;i<ysiz;i++)
		{
			mptr[i] = data + i*xsiz;
		}

		return mptr;
	}
	catch( std::bad_alloc &ba )
	{
		std::string mess = "AllocMatrix : メモリの確保に失敗！\r\n";
		mess += ba.what();
		Message(mess.c_str());
		return NULL;
	}
}

// メモリ解放
template <typename T> void Free(T **mptr)
{
	if( mptr )
	{
		delete [] mptr[0];
		delete [] mptr;
	}
}

template <typename T> bool LoadMatrix( T **array , int xsize, int ysize , FILE *fp )
{
	T data;

	for(int m=0 ; m<ysize ; m++)
	{
		for(int n=0 ; n<xsize ; n++)
		{
			if( fscanfVar( &data , fp ) != 1 )
			{
				return false;
			}
			fscanf( fp , "	" );

			array[m][n]=data;
		}
		
		fscanf( fp , "\n");
	}
	return true;
}


template <typename T> T **LoadMatrix(int &xsize, int &ysize , const char *filename )
{
	FILE *fp = NULL;
	T **array = NULL;
	char mess[1024];

	fp = fopen( filename , "r");

	if( fp==NULL )
	{
		sprintf( mess, "%s がありません．" , filename );
		Message( mess );

	}


	//行数取得
	ysize = CountLine(fp);
	xsize = CountCols(fp);

	array = (T **)AllocMatrix<T>( ysize , xsize );

	if( !LoadMatrix( array , xsize, ysize , fp ) )
	{
		if( !fp ) fclose( fp );
		if( !array ) Free( array );

		sprintf( mess, "%s の書式が不正です．" , filename );
		Message( mess );

		return NULL;
	}

	fclose(fp);

	return array;
}