/** applet No. 1044
 *
 * harmonic oscillator - path-integral quantum Monte-Carlo method 1D
 * - multi thread  PIMC :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on October 22 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.10.22  created
 * ver 0.0.2   2006.12.04  improved code
 * ver 0.0.3   2007.05.26  improved code
 *
 */

import java.awt.*;
import java.awt.event.*;
import java.applet.*;

public class harmonicsMtPIMC1D extends Applet
  implements ItemListener, ActionListener, AdjustmentListener, Runnable {

  // ------------------------------------ preset field -----------

  Thread th = null; // for run()-paint() thread
  PathIntegralMC pimc = null;

  Choice ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_nPath, sc_kConst;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int ithPath = 0;
  int dispMode = 2;

  int sleepTime = 30;
  int thCount = 0;


  // ------------------------------ applet main thread -----------

  public void init() {

    resize(630,320);
    setBackground(Color.white);

    dim = getSize();
    imgOff = createImage(dim.width,dim.height);
    gOff = imgOff.getGraphics();

    ch_view = new Choice();
    ch_view.add("density");
    ch_view.add("wave");
    ch_view.add("path");
    ch_view.addItemListener(this);
    ch_view.select("path");

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 

    bt_step= new Button("step");
    bt_step.addActionListener(this); 

    sc_nPath = new Scrollbar(Scrollbar.HORIZONTAL,100,10,1,210);
    sc_nPath.addAdjustmentListener(this);

    sc_kConst = new Scrollbar(Scrollbar.HORIZONTAL,50,10,50,210);
    sc_kConst.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(bt_step);
    pnl.add(sc_nPath);
    pnl.add(sc_kConst);
    pnl.add(ch_view);
    add(pnl,"North");

  }

  public void start() {
    if (pimc == null) {
      pimc = new PathIntegralMC();
      pimc.start();
    }
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

  public void stop() {
    if (pimc == null) pimc = null;
    if (th == null) th = null;
  }


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      pimc.reset();
      thCount = 0;
    } else if (ev.getSource() == bt_startStop){
      if (pimc.getStartSW()==0) {
        pimc.setStartSW(1);
      } else {
        pimc.setStartSW(0);
      }
    } else if (ev.getSource() == bt_step){
      if (pimc.getStartSW()==0) pimc.setStepSW();
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_nPath) { 
      pimc.setnMax( sc_nPath.getValue() );
    } else if (ev.getSource() == sc_kConst) { 
      pimc.setkConstant( 0.01*(double)(sc_kConst.getValue()) );
    }
  }



  // ========================= run() - paint() loop ==============

  public void run() {
    while (th != null) {
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
      thCount += 1;
      offPaint();
      repaint();
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }


  // ---------------------------------------- offPaint -----------

  void offPaint() {
    int gbx,gby;
    double sc,dispScale,viewScale;

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,dim.width,dim.height);

    gbx = 30;
    gby = 60;
    dispScale = (250.0/(10.0e-9));
    viewScale = 1.0;
    sc = viewScale*dispScale;

    if (dispMode==0) {
      densityPlot(gbx,gby, 5000.0);
      potentialPlot(gbx, gby);

    } else if (dispMode==1) {
      wavePlot(gbx,gby,500.0);
      potentialPlot(gbx, gby);

    } else if (dispMode==2) {
      wavePlot(gbx,gby,500.0);
      potentialPlot(gbx, gby);
      pathPlot(gbx, gby, ithPath);

    } else if (dispMode==3) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("PIMC="+pimc.getCount()+" ",10,40);
    if (pimc.getStartSW()==0) {
      gOff.drawString("stopped",630/6*1+10,40);
    } else {
      gOff.drawString("started",630/6*1+10,40);
    }
    gOff.drawString("Npath="+pimc.getnMax()+"",630/6*3+10,40);

    gOff.drawString("k="+pimc.getkConstant()+"",630/6*4+10,40);
    gOff.drawString("V(x) = k x^2",630/6*4+10,60);
    gOff.drawString("view",630/6*5+10,40);

    gOff.drawString("norm ="+(float)pimc.getNorm()+"",400,100);
    gOff.drawString("<psi|H|psi> ="+(float)pimc.getEnergy()+"",400,120);
    gOff.drawString("virial energy ="+(float)pimc.getMeanVirialEnergy()+"",400,140);

    gOff.drawString("loop count",400,200);
    gOff.drawString("thread disp ="+thCount+"",400,220);
    gOff.drawString("thread pimc ="+pimc.getCount()+"",400,240);
    gOff.drawString("disp sleep = 30 ms",400,260);

  }


  // ------------------------------------ plot methods -----------

  void densityPlot(int gbx, int gby, double mag) {
    int ix, iy, dmax;

    gOff.setColor(Color.getHSBColor(0.66f,0.3f,0.9f ));
    dmax = pimc.getdMax();
    for(ix=0; ix<dmax; ix++) {
      iy = (int)(mag*pimc.getMeanDensity(ix));
      gOff.fillRect(gbx+2*ix,gby+200-iy,2,iy);
    }

  }

  void wavePlot(int gbx, int gby, double mag) {
    int ix, iy, dmax;

    gOff.setColor(Color.getHSBColor(0.33f,0.3f,0.9f ));
    dmax = pimc.getdMax();
    for(ix=0; ix<dmax; ix++) {
      iy = (int)(mag*pimc.getwaveFunction(ix));
      gOff.fillRect(gbx+2*ix,gby+200-iy,2,iy);
    }
  }

  void potentialPlot(int gbx, int gby) {
    int ix, iy, iy2,dmax;

    gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.9f ));
    dmax = pimc.getdMax();
    for(ix=0; ix<dmax-1; ix++) {
      iy = (int)(20.0*pimc.getPotential(ix));
      iy2 = (int)(20.0*pimc.getPotential(ix+1));
      gOff.drawLine(gbx+2*ix,gby+200-iy,gbx+2*ix+2,gby+200-iy2);
    }
  }

  void pathPlot(int gbx, int gby, int ipath) {
    int j, ix, ix2,dmax, mmax;
    double x,x2;

    gOff.setColor(Color.getHSBColor(0.15f,0.3f,0.8f ));
    dmax = pimc.getdMax();
    mmax = pimc.getmMax();

    for(j=0; j<mmax-1; j++) {
      x = pimc.getPath(ipath,j);
      ix = (int)(2.0*(x-pimc.xMin)/pimc.dx);
      x2 = pimc.getPath(ipath,j+1);
      ix2 = (int)(2.0*(x2-pimc.xMin)/pimc.dx);
      gOff.setColor(Color.getHSBColor((float)(j*1.0/mmax),0.3f,0.8f ));
      gOff.drawLine(gbx+ix,gby+2*j,gbx+ix2,gby+2*j+2);
      gOff.setColor(Color.getHSBColor(0.15f,0.8f,0.8f ));
      gOff.drawOval(gbx+ix-3,gby+2*j-3,6,6);
    }
  }
}



// ============================ path-integral MC class ===========

class PathIntegralMC extends Thread {

  public static final double xMin = -4.0;
  public static final double xMax = 4.0;
  public static final double dx = 0.05;

  private int dMax = (int)((xMax-xMin)/dx+1);

  private int nMax = 100;
  private int mMax = 100;
  private double tau = 150; // imaginal time
  private double dtau = tau/mMax;
  private double dt2 = 2.0*dtau*dtau;
  private double temperature = 1/tau;
  private double aWidth = (xMax-xMin)/10.0;

  private double kConstant = 0.5;

  private int NNN = 1000;
  private int MMM = 100;
  private double path[][] = new double[NNN][MMM];

  private int nDensity[] = new int[dMax];
  private double density[] = new double[dMax];
  private double meanDensity[] = new double[dMax];
  private double waveFunction[] = new double[dMax];
  private double meanVirialEnergy = 0.0;


  private int resetSW = 0; // 1-reset
  private int startSW = 1; // 0-stop 1-start
  private int stepSW = 0;
  private int sleepTime = 1;
  private int waitSleepTime = 100;
  private int count = 0;


  // ------------------------------- class constructor -----------

  PathIntegralMC() {
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      count += 1;
      timeEvolution();
      try {
        this.sleep(sleepTime);
      } catch (InterruptedException e) { }
    }
  }


  // --------------------------- set initial condition -----------

  void setInitialCondition() {
    int i,j;

    count = 0;

    for (i=0; i<nMax; i++) {
      for (j=0; j<mMax; j++) {
        path[i][j] = 2.0*(Math.random()-0.5);
      }
    }

    for (i=0; i<dMax; i++) {
        meanDensity[i] = 0.0;
    }
    meanVirialEnergy = 0.0;
  }


  // ---------------------------------- time evolution -----------

  void timeEvolution() {
    int i;

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }

    if (startSW==1) {
        timeStep();
    } else {
      if (stepSW==1) {
        timeStep();
        stepSW = 0;
      }
      try {
        this.sleep(waitSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void timeStep() {
    int i,j;

    for (i=0; i<nMax; i++) {
      for (j=0; j<mMax; j++) {
        samplingPath(path[i]);
      }
    }
    setDensity();
    setMeanVirialEnergy();
  }

  // path-integral sampling

  private int samplingPath(double xx[]) {
    int j,jp,jm;
    double xtrial,dEnergy;

    j = (int)(mMax*Math.random()); if (j==mMax) j = mMax-1;
    jp = (j+1)%mMax; jm = (j-1+mMax)%mMax;

    do {
      xtrial = xx[j] + aWidth*(2.0*Math.random()-1.0);
    } while ( xtrial<xMin || xtrial>xMax );
    dEnergy = ( (xx[jp]-xtrial)*(xx[jp]-xtrial)+(xtrial-xx[jm])*(xtrial-xx[jm])
      -(xx[jp]-xx[j])*(xx[jp]-xx[j]) -(xx[j]-xx[jm])*(xx[j]-xx[jm]) )/dt2
	+ potential(xtrial) - potential(xx[j]);

    if ( Math.exp(-dEnergy)>Math.random() ) {
      xx[j] = xtrial;
      return j;
    }
    return -1;
  }

  private double potential(double x) {
    return ( kConstant*x*x );
  }

  // set density

  private void setDensity() {
    int i,j,ix,nm;
    double x,ei;

    nm = nMax*mMax;
    ei = 1.0/(count+1.0); if (count>10000) ei=0.0001;

    for (ix=0; ix<dMax; ix++) {
      nDensity[ix] = 0;
    }

    for (i=0; i<nMax; i++) {
      for (j=0; j<mMax; j++) {
        x = path[i][j]-xMin;
        ix = (int)(x/dx);
        nDensity[ix] += 1;
      }
    }

    for (ix=0; ix<dMax; ix++) {
      density[ix] = (1.0*nDensity[ix])/nm;
      meanDensity[ix] = (1.0-ei)*meanDensity[ix]+ei*density[ix];
      waveFunction[ix] = Math.sqrt(meanDensity[ix]);
    }
  }

  //

  private void setMeanVirialEnergy() {
    double ei;

    ei = 1.0/(count+1.0); if (count>10000) ei=0.0001;
    meanVirialEnergy = (1.0-ei)*meanVirialEnergy + ei*virialEnergy();
  }

  private double virialEnergy() {
    int i;
    double s;

    s = 0.0;
    for (i=0; i<nMax; i++) {
      s += virialEnergy(i);
    }
    return s/nMax;
  }

  private double virialEnergy(int i) {
    int j;
    double s,x,h;

    h = 0.01;
    s = 0.0;
    for (j=0; j<mMax; j++) {
      x = path[i][j];
      s += potential(x) + 0.5*x*(potential(x+h)-potential(x-h)/(2.0*h));
    }
    return s/mMax;
  }


  // ----------------------------------------- utility -----------

  private double norm() {
    int ix;
    double s;

    s = 0.0;
    for (ix=0; ix<dMax; ix++) {
      s += waveFunction[ix]*waveFunction[ix];
    }
    return s;
  }

  private double energy() {
    int ix;
    double s;

    s = 0.0;
    for (ix=1; ix<dMax-1; ix++) {
      s += -0.5*waveFunction[ix]*(waveFunction[ix+1]+waveFunction[ix-1]-2.0*waveFunction[ix])/(dx*dx);
      s += waveFunction[ix]*waveFunction[ix]*potential(ix*dx+xMin);
    }
    return s;
  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; }
  public void setnMax(int nm) { nMax = nm; resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public int getCount() { return count; }

  public void setkConstant(double kC) { kConstant = kC; }
  public double getkConstant() { return kConstant; }

  // get data
  public int getnMax() { return nMax; }
  public int getmMax() { return mMax; }
  public double getPath(int i, int j) { return path[i][j]; } 

  public int getdMax() { return dMax; }
  public double getDensity(int ix) { return density[ix]; }
  public double getMeanDensity(int ix) { return meanDensity[ix]; }
  public double getwaveFunction(int ix) { return waveFunction[ix]; }
  public double getNorm() { return norm(); }
  public double getEnergy() { return energy(); }
  public double getMeanVirialEnergy() { return meanVirialEnergy; }
  public double getPotential(int ix) { return potential(ix*dx+xMin); }

}


/** path-integral quantum Monte-Carlo method
 *
 * - classical Metropolis random sampling in N+1 dimension space
 *
 * time evolution operator in imaginary time 
 *  exp(-tau H)
 *
 * partition function:
 * Tr(exp(-bata H)) = Integrate(<R|exp(-bata H)|R>,dR)
 *
 * Integrate(<R0|exp(-bata H)|R0>,dR0) 
 *  = Integrate(<R0|exp(-bata H)|R1><R1|exp(-bata H)|R2><R2..
 *      ..R_M-1><R_M-1|exp(-bata H)|R0>,dR0...dR_M-1 )
 *
 * time evolution in short time:
 * Integrate(<R|exp(-tau H)|R'>,dR0)
 *  ~ normalized exp(-dtau V(R)) exp(-(R-R')^2/2dtau)
 *
 *
 * Integrate(<R0|exp(-tau H)|R0>,dR0)
 *  ~ normalized Integrate( exp(-dtau*sum(Pm,{m,0,M-1})))
 *  ( where Pm = {(R_m+1-Rm)/dtau }^2 + V(Rm) )
 *
 *
 * method
 *
 * path - {x[i] ,i=0,..,j,..,M-1}
 *   - sampling path (Metropolis Algorithm)
 *
 * (1) trial move at j
 *  select j at random:
 *   trial move x[j] + xtrial, xtrial <-- a*random[-1,1]
 *
 * (2) energy change : dEnergy
 *  dEnergy = { (x[j+1]-xtrial)^2 + (xtrial-x[j-1])^2 
 *   - (x[j+1]-x[j])^2 - (x[j]-x[j-1])^2 }/(2dtau^2) 
 *   + {V(xtrial)-V(x[j]) }
 *
 * (3) move or not
 *  if exp(-dEnergy)>random[0,1] 
 *    then move x[j] <-- xtrial
 *
 * goto (1)
 *
 */