artmap.cpp

ARTMAP的C语言程序

    
    for ( lambda=0; lambda<C; lambda++ )
      Sum += y[lambda];
    
    for ( j=0; j<C; j++)
      Y[j] = y[j]/Sum;
  }
  
  // F3 -> F1 signal
  for(i=0; i<2*M; i++)
    {
      Sum = 0;
      for(j=0; j<C; j++ )
	Sum += max(Y[j]-tau_ji[j][i], 0);
      sigma_i[i] = Sum;
    }            
  
 Step_8: //Resonance
  for( j=0; j<C; j++ )
    {
      for( i=0; i< 2*M; i++ )
	{
	  tau_ij[i][j] = tau_ij[i][j] + beta * max(y[j]-tau_ij[i][j]-A[i], 0);
	  aux = beta * max(sigma_i[i]-A[i], 0) * max( Y[j] - tau_ji[j][i], 0 );
	  
	  if( sigma_i[i] != 0.0 ) aux = aux/sigma_i[i];
	    
	  tau_ji[j][i] = tau_ji[j][i] + aux;
        }
      c[j] = c[j] + y[j];
    }
  Delta_len = 0; //reset node recovery
  rho = rho_a_bar; //ARTa vigilance recovery
  
 Step_9: // Next iteration
  if(n==TRAIN_N-1)
    {
      //Training_output
      printf("Epoch: %3d     Commited F2 nodes: %3d\n",epochs+1, C);
      if(epochs==EPOCHS-1)
	return; // End of training
      epochs++;
      n=-1;
    }
  n++;
  
  // Copy next input pattern and corresponding output category
  for( i=0; i<M; i++ )
    {
      A[i] = input[n][i];
      A[i+M] = 1-A[i];
    }
  K = output[n];
  dist_mode = 1; //revert to distributed mode
  goto Step_2;
}

void test()
{
  printf("Testing\n");
  cnum = 0;
  
 Test_Step_1: //First iteration
  for(test_n = 0; test_n < TEST_N; test_n++)
    {
      // Copy next input pattern and corresponding output category
      for( i=0; i<M; i++ )
        {
	  A[i] = te_input[test_n][i];
	  A[i+M] = 1-A[i];
	}
      K = te_output[test_n];
      
    Test_Step_2:  // Reset
      bu_match();
      
      if ( DO_TEST_ICG ) {
	Lambda_len = 0;
	for( j=0; j<C; j++ )
	  if( T[j] >= T_u )
	    Lambda[Lambda_len++] = j;
	
	// (a) If the network is in distributed mode: F2 nodes are activated
	//     according to the increased gradient CAM rule.
	Lambda_pp_len = 0;
	for(j=0; j<Lambda_len; j++)
	  if( M - T[Lambda[j]] == 0 )
	    Lambda_pp[Lambda_pp_len++] = Lambda[j];
	
	for( j=0; j<C; j++ )
	  y[j]=0;
	
	// (i) If M-Tj>0 for all j belonging to Lambda...
	if( Lambda_pp_len == 0 )
	  for( j=0; j<Lambda_len; j++)
	    {
	      Sum = 0;
	      for(lambda=0;lambda<Lambda_len;lambda++)
                {
		  if(Lambda[lambda]==Lambda[j])   continue;
		  Sum += pow((M-T[Lambda[j]]) / 
			     (M-T[Lambda[lambda]]),p);
                }
	      y[Lambda[j]] = 1./(1.+Sum);
	    }
	// (ii) Point box case:
	else
	  for( j=0; j< Lambda_pp_len; j++ )
	    y[Lambda_pp[j]] = 1./Lambda_pp_len;
	
	// F3 activation
	if ( DO_TEST_IC )
	  {
            double Sum_clyl = 0;
            for( lambda=0; lambda<C; lambda++ )
	      
	      Sum_clyl += c[lambda] * y[lambda];
            for( j=0; j<C; j++ )
	      Y[j] = c[j]*y[j] / Sum_clyl;
	  }
	else {
	  double Sum_clyl = 0;
	
	  for ( lambda=0; lambda<C; lambda++ )
	    Sum_clyl += y[lambda];

	  for ( j=0; j<C; j++)
	    Y[j] = y[j]/Sum_clyl;
	}
      }

      if ( DO_TEST_SCG ) {
	Sum = 0;
	for ( i = 0; i < C; i++)
	  Sum += pow(T[i], p);

	for( i = 0; i < C; i++)
	    T[i] = pow(T[i],p)/Sum;

	// F2 activation
	Sum = 0;
	for ( i = 0; i < C; i++)
	  Sum += T[i];

	for ( i = 0; i < C; i++)
	  y[i] = T[i]/Sum;
	
	// F3 activation
	if ( DO_TEST_IC ) {
	  double Sum_clyl = 0;
	  for( lambda=0; lambda<C; lambda++ )
	    Sum_clyl += c[lambda] * y[lambda];
	  
	  for( j=0; j<C; j++ )
	    Y[j] = c[j]*y[j] / Sum_clyl;
	}
	else {
	  double Sum_clyl = 0;
	  
	  for ( lambda=0; lambda<C; lambda++ )
	    Sum_clyl += y[lambda];
	  
	  for ( j=0; j<C; j++)
	    Y[j] = y[j]/Sum_clyl;
	}
      }

      if ( DO_TEST_WTA ) {
	// (b) If the network is in WTA mode: Only one F2 node, with j=J, is activated
	// (i) If there is a commited node:
	if( Lambda_len > 0 )
	  {
	    J = Lambda[0];
	    for( j=1; j<Lambda_len; j++ )
	      if( T[Lambda[j]] > T[J]) 
		J=Lambda[j];
	  }
	else    // No "commited" nodes, but pick the best you can
	  {
	    double Tmax = -1.0;
	    for (j = 1; j < C; j++) {
	      if ( T[j] > Tmax ) {
		Tmax = T[j];
		J = j;
	      }
	    }
	  }
	
	// (ii) F2 and F3 activation
	for( j=0; j<C; j++ )
	  y[j] = Y[j] = 0.;
	y[J] = Y[J] = 1.;
      }

Test_Step_3: // Prediction
        for( k=0; k<L; k++ )
	  sigma_k[k] = 0;
        for( j=0; j<C; j++ )
	  sigma_k[kappa[j]] += Y[j];
	
        K_prime = 0;
        for( k=1; k<L; k++ )
	  if( sigma_k[k] > sigma_k[K_prime] )
	    K_prime = k;

Test_Step_4: // Evaluation
        if( K_prime == K )
	  {
            cnum++;
	  }
        else
	  //Testing_output
	  printf("Test pat #%4d: INcorrect prediction\n", test_n+1);
    }   
  printf("Number of correctly classified test patterns: %d\n",cnum);
}

int is_in_Delta( int j, int *Delta, int Delta_len)
{
  int in_Delta = 0;
  
  for(int  in_Delta_index = 0; in_Delta_index < Delta_len; in_Delta_index++ )
    if( j == Delta[in_Delta_index] )
      {
            in_Delta = 1;
            break;
      }
  return in_Delta;
}

void getInputOutput(char* train_input, char* train_output, char* test_input, char* test_output) {
    // read training and testing data from files
    FILE* input_f     = fopen(train_input,"r");
    FILE* output_f    = fopen(train_output,"r");
    FILE* te_input_f  = fopen(test_input,"r");
    FILE* te_output_f = fopen(test_output,"r");
    int i, j;

    if( ! (input_f && output_f && te_input_f && te_output_f) )
    {
      printf("Error: Unable to open one of the data files!\n");
      exit(0);
    }
    
    for( i=0; i<TRAIN_N; i++ )
      for(int j=0; j<M; j++ )
	fscanf(input_f,"%lf",&(input[i][j]));
    
    for( i=0; i<TRAIN_N; i++ )
      {
        fscanf(output_f,"%d",&output[i]);
        output[i] = output[i] - 1;
      }
    
    for( i=0; i<TEST_N; i++ )
      for(int j=0; j<M; j++ )
	fscanf(te_input_f,"%lf",&(te_input[i][j]));
    
    for( i=0; i<TEST_N; i++ )
      {   
        fscanf(te_output_f,"%d",&(te_output[i]));  
        te_output[i] = te_output[i] - 1;
      }    
    
    fclose(input_f);
    fclose(output_f);
    fclose(te_input_f);
    fclose(te_output_f);
}

void forceHypercube() {
  double max[M], min[M];
  for ( int i = 0; i < M; i++ ) {
    max[i] = input[0][i];
    min[i] = input[0][i];
  }

  for ( int i = 0; i < TRAIN_N; i++ ) {
    for ( int j = 0; j < M; j++) {
      if (input[i][j] > max[j] )
	max[j] = input[i][j];

      if (input[i][j] < min[j] )
	min[j] = input[i][j];
    }
  }

  for ( int i = 0; i < TRAIN_N; i++ ) {
    for ( int j = 0; j < M; j++) {
      input[i][j] = (input[i][j] -  min[j]) / (max[j]-min[j]);
    }
  }

  for ( int i = 0; i < M; i++ ) {
    max[i] = te_input[0][i];
    min[i] = te_input[0][i];
  }

  for ( int i = 0; i < TEST_N; i++ ) {
    for ( int j = 0; j < M; j++) {
      if (te_input[i][j] > max[j] )
	max[j] = te_input[i][j];

      if (te_input[i][j] < min[j] )
	min[j] = te_input[i][j];
    }
  }

  for ( int i = 0; i < TEST_N; i++ ) {
    for ( int j = 0; j < M; j++) {
      te_input[i][j] = (te_input[i][j] -  min[j]) / (max[j]-min[j]);
    }
  }

}

void checkInputOutput() {
  int quit = 0;

  for ( int i = 0; i < TRAIN_N && !quit; i++ ) {
    for ( int j = 0; j < M && !quit; j++) {
      if (input[i][j] > 1.0 || input[i][j] < 0.0)
	quit = 1;
    }
  }

  if ( quit == 1 ) {
    printf( "Training input is not in unit hypercube!\n" );
    exit(0);
  }

  quit = 0;

  for ( int i = 0; i < TEST_N && !quit; i++ ) {
    for ( int j = 0; j < M && !quit; j++) {
      if (te_input[i][j] > 1.0 || te_input[i][j] < 0.0)
	quit = 1;
    }
  }

  if ( quit == 1 ) {
    printf( "Testing input is not in unit hypercube!\n" );
    exit(0);
  }
}

void bu_match() {
  for( j=0; j<C; j++ )
    {
      if (is_in_Delta( j, Delta, Delta_len) )
	T[j] = 0;
      else
        {
	  if ( DO_CBD ) {
            // Choice by difference function
            S[j] = 0; Theta[j] = 0;
            for( i=0; i<2*M; i++ )
	      {
                S[j] += min(A[i], 1-tau_ij[i][j]);
                Theta[j] += tau_ij[i][j];
	      }
            T[j] = S[j]+(1-alpha)*(Theta[j] - M);
	  }
	  
	  if ( DO_WEBER ) {
            // Weber's law choice function
            T[j]=Sum=0;
            for( i=0; i<2*M; i++ )
	      {
                T[j] += min(A[i], 1-tau_ij[i][j]);
                Sum += 1-tau_ij[i][j];
	      }
            T[j] = T[j] / (alpha + Sum );
	    
	    T_u = M/(alpha + 2.0*M);
	  }

        }
    }
}

void printArr( char* str, double* arr, int len) {
  int i;

  fprintf( stdout, "**\t%s = ", str );
  for (i = 0; i < len-1; i++) {
    fprintf( stdout, "%5.2f, ", arr[i] );
  }
  fprintf( stdout, "%5.2f\n", arr[i] );
}

void printArra( char* str, int* arr, int len) {
  int i;

  fprintf( stdout, "**\t%s = ", str );
  for (i = 0; i < len-1; i++) {
    fprintf( stdout, "%d, ", arr[i] );
  }
  fprintf( stdout, "%d\n", arr[i] );
}

void printArr2( char* str, double arr[F0_SIZE][F2_SIZE], int len1, int len2) {
  int i, j;

  for (i = 0; i < len1; i++) {
    fprintf( stdout, "**\t%s[%d] = ", str, (i+1) );
    for (j = 0; j < len2-1; j++) {
      fprintf( stdout, "%5.2f, ", arr[i][j] );
    }
    fprintf( stdout, "%5.2f\n", arr[i][j] );
  }
}

void printArr2a( char* str, double arr[F2_SIZE][F0_SIZE], int len1, int len2) {
  int i, j;

  for (i = 0; i < len1; i++) {
    fprintf( stdout, "**\t%s[%d] = ", str, (i+1) );
    for (j = 0; j < len2-1; j++) {
      fprintf( stdout, "%5.2f, ", arr[i][j] );
    }
    fprintf( stdout, "%5.2f\n", arr[i][j] );
  }
}