#!/usr/bin/env --split-string=${JDK_HOME}/bin/java --source ${JDK_VERSION}
package scratch;

import static java.lang.Double.parseDouble;
import static java.lang.Math.pow;
import static java.lang.System.err;
import static java.lang.System.exit;
import static java.lang.System.out;


/** A calculator of the likelihood (P) of an event of decaying probability.  The calculation is a sum
  * over a large number of trials (t = 1, 2, 3, …) using two parameters: the initial per-trial
  * probability 0 ≤ p ≤ 1 of the event; and the per-trial decay rate 0 ≤ δ ≤ 1 of the probability
  * (applied after each trial), where δ = 1 yields a constant probability. *//*
  *
  *   https://math.stackexchange.com/questions/4692693/bernoulli-series-with-decaying-probability-p
  */
class P {


    public static void main( final String[] arguments ) {
        exitOnDemand( arguments );
        final int n = arguments.length;
        if( n != 2 ) {
            err.println( commandName + ": Expecting 2 arguments, found " + n );
            exitWithUsage( 1 ); }
        new P().run( /*p*/parseDouble(arguments[0]), /*δ*/parseDouble(arguments[1]) ); }



////  P r i v a t e  ////////////////////////////////////////////////////////////////////////////////////


    static final String commandName = "P";



    static void exitOnDemand( final String[] arguments ) {
        for( final String arg: arguments ) {
            if( arg.equals("-?") || arg.equals("-help") ) exitWithUsage( 0 ); }}



    static void exitWithUsage( final int status ) {
        out.println( commandName +
          ": Likelihood of an event of initial probability p, decaying at rate δ" );
        out.println( "Usage: " + commandName + " p δ" );
        out.println( "       " + commandName + " -help | -?" );
        exit( status ); }



    /** Cumulative likelihood of occurence.
      */
    double pCumulative;



    /** Per-trial probability at last trial.
      */
    double pLast;



    /** Cumulative likelihood of failure.
      */
    double qCumulative = 1d;



    /** @param p Per-trial probability, initial value.
      * @param δ Per-trial decay rate of `p`.
      */
    void run( double p, final double δ ) {
        testParameter( p );
        testParameter( δ );
        final double δInverse = 1d - δ;
        final double p0 = p;
        do {
            qCumulative *= 1d - pLast;
            final double pJustNow/*earlier failure *and* occurence on the present trial*/ =
              qCumulative/*likelihood of earlier failure*/ * p/*per-trial probability*/;
            pCumulative += pJustNow;
            if( t == tReport ) {
                out.println( "t " + t + "     \t" + pCumulative );
                if( tReport < 20 ) tReport += 1;
                else if (tReport == 20 ) tReport = trialReportInterval;
                else tReport += trialReportInterval; }
            pLast = p;
            p *= δInverse; }
            while( ++t <= trialsMaximum );
        out.println();
        out.println( "  % p0 = " + p0 * 100d + "                (p initial)" );
        out.println( "       = " + (1d - pow(1d-p0,100)) * 100d + " / century" );
        out.println( "   % δ = " + δ * 100d );
        out.println( "       = " + (1d - pow(δInverse,100)) * 100d + " / century" );
        out.println( "   % P = " + pCumulative * 100d + "  (extinction)" );
        out.println( "   % Q = " + (1d - pCumulative) * 100d + "  (survival)" ); }



    /** Trial number.
      */
    int t = 1;



    static void testParameter( final double d ) {
        if( 0d > d || d > 1d ) throw new IllegalArgumentException( "Parameter outside range 0 to 1" ); }



    /** Trial number for the next report.
      */
    int tReport = t;



    static final int trialReportInterval = 10000;



    static final int trialsMaximum = trialReportInterval * 12; }