printing results as txt

RUN pip3 install scipy numpy matplotlib lmfit

RUN pip3 install --upgrade tensorflow

RUN pip3 install tensorrt

RUN pip3 install numdifftools

COPY /src /src

<tool id="mlsans_fit_Yukawa_sq" name="fit_IQ"  profile="22.04" version="0.1.1">

  <description></description>

    <requirements>

        <container type="docker">code.ornl.gov:4567/ndip/tool-sources/ml-assisted-sans-data-analysis/fit_yukawa_sq:99dee52e</container>

        <container type="docker">code.ornl.gov:4567/ndip/tool-sources/ml-assisted-sans-data-analysis/fit_yukawa_sq:84114dfe</container>

    </requirements>

    <command detect_errors="exit_code"><![CDATA[

        python3 /src/curvefit.py -i $input $C $I_inc $sigma $d_sigma $phi $kappa $A 

        python3 /src/curvefit.py -i $input $C $I_inc $sigma $d_sigma $phi $kappa $A || echo python finished

    ]]></command>

    <inputs>

        <param name="input" type="data" optional="false" label="ASCII data"/>

        <param name="input" type="data" value="/src/example_IQ.csv" optional="false" label="ASCII data"/>

        <param name="C" type="text" value="5000, 1000, 8000" label="C" optional="false" help="Contrast"/>

        <param name="I_inc" type="text" value="0.2, 0.1, 0.5" label="I_inc" optional="false" help="Incoherent background"/>

        <param name="sigma" type="text" value="1, 0.9, 1.1" label="σ" optional="false" help="Particle diameter"/>

    </inputs>

    <outputs>

        <data name="output" format="png" label="LogLog Plot">

        <data name="output" format="txt" label="fitting results">

        </data>

    </outputs>

    <help><![CDATA[

%% Cell type:code id: tags:

``` python

import numpy as np

import matplotlib.pyplot as plt

import scipy.interpolate as interp

from scipy.signal import savgol_filter

from lmfit import Minimizer, Parameters, create_params, report_fit, Model

import os

```

%% Cell type:code id: tags:

``` python

from SQ_NN import SQ_NN, SQ_NN_tf, Decoder_aug

QD = (np.arange(100)+1)*0.2

```

%% Cell type:code id: tags:

``` python

def hardsphere(q, sigma=1):

    R = sigma / 2

    P = (3 * (np.sin(q * R) - q * R * np.cos(q * R)) / (q * R) ** 3) ** 2

    return P

def fuzzysphere(q, sigma=1, sigma_f=0.1):

    R = sigma / 2

    P = (3 * (np.sin(q * R) - q * R * np.cos(q * R)) / (q * R) ** 3) ** 2 * np.exp(-(sigma_f * sigma * q) ** 2 / 2)

    return P

def log_normal_pdf(mu, sigma, x):

    return np.exp(-(np.log(x) - mu) ** 2 / 2 / sigma ** 2) / x / sigma

def P_HS_eff(q, sigma=1, d_sigma=0.05, return_f=False):

    '''

    sigma: average particle size

    d_sigma: polydispersity

    return_f: toggle whether to return the particle size distribution

    '''

    # List of particle diameter

    n_sample = 101

    sigma_list = (1 + np.linspace(-5, 5, n_sample) * d_sigma) * sigma

    sigma_list = sigma_list[sigma_list > 0]

    # Size distribution

    f_sigma = log_normal_pdf(0, d_sigma, sigma_list / sigma)

    p_sigma = f_sigma * (sigma_list / sigma) ** 6

    # Calculate effective P(Q)

    P_eff = np.zeros_like(q)

    for i in range(len(sigma_list)):

        P_i = hardsphere(q, sigma_list[i]) * p_sigma[i]

        P_eff = P_eff + P_i

    P_eff = P_eff / np.sum(p_sigma)

    if return_f:

        return P_eff, sigma_list, f_sigma

    else:

        return P_eff

def IQ_th(params, Q):

    v = params.valuesdict()

    fp = [v['phi'],v['kappa'],v['A'],1]

    # structure factor

    S = SQ_NN(fp[0:3])

    S = savgol_filter(S,7,2)

    # form factor

    P = P_HS_eff(Q,sigma=v['sigma'], d_sigma=v['d_sigma'])

    I = v['C']*P*S + v['I_inc']

    return I

```

%% Cell type:code id: tags:

``` python

## Generate example input

# initialize parameters

params = Parameters()

params.add('C', value=4100, min=4000, max=5000)

params.add('I_inc', value=0.5, min=0.55, max=2)

params.add('I_inc', value=0.5, min=0.1, max=2)

params.add('sigma', value=0.975, min=0.9, max=1.1)

params.add('d_sigma', value=0.07, min=0.05, max=0.1)

params.add('phi', value=0.08, min=0.05, max=0.5)

params.add('kappa', value=0.15, min=0.1, max=0.18)

params.add('A', value=10, min=0.5, max=20)

IQ_example = IQ_th(params, QD)

import csv

path = './'

filename = 'example_IQ.csv'

with open(path+filename, 'w', newline='') as f:

    writer = csv.writer(f)

    for (QD_i, IQ_i) in zip(QD,IQ_example):

        writer.writerow([QD_i, IQ_i])

```

%% Cell type:code id: tags:

``` python

data = np.genfromtxt(path+filename,delimiter=',')

QD_exp = data[:,0]

IQ_exp = data[:,1]

fig = plt.figure(figsize=(6, 6))

ax = plt.subplot(1, 1, 1)

ax.plot(QD_plot,IQ_plot,'.b',label='SANS data')

ax.set_yscale('log')

ax.set_xscale('log')

# ax.set_xlim([2,40])

# ax.set_ylim([0.2,4000])

ax.set_xlabel(r'$QD$',fontsize=20)

ax.set_ylabel(r'$I(QD)$',fontsize=20)

ax.tick_params(direction='in', axis='both', which='both', labelsize=18, pad=10)

ax.legend(fontsize=16,frameon=False)

plt.show()

```

%% Output

%% Cell type:code id: tags:

``` python

## Test fitting code

# initialize parameters

params = Parameters()

params.add('C', value=5000, min=4000, max=6000)

params.add('I_inc', value=0.4, min=0.1, max=2)

params.add('sigma', value=1, min=0.9, max=1.1)

params.add('d_sigma', value=0.07, min=0.05, max=0.1)

params.add('phi', value=0.2, min=0.05, max=0.5)

params.add('kappa', value=0.2, min=0.1, max=0.18)

params.add('A', value=10, min=0.5, max=20)

def lmbda(params, Q, IQ_exp, index_Q):

    IQ = IQ_th(params, Q)

    minimizer_target = lambda x, y, z: (np.log(x/y))/np.log(1+z/y)

    # minimizer_target = lambda x, y, z: (x-y)**2

    return minimizer_target(IQ[index_Q],IQ_exp[index_Q],np.ones_like(IQ_exp))

# do fit, here with the nelder algorithm

minner = Minimizer(lmbda, params, fcn_args=(QD_exp, IQ_exp, np.arange(len(QD_exp))))

result = minner.minimize('nelder')

# write error report

report_fit(result)

```

%% Output

    [[Fit Statistics]]

        # fitting method   = Nelder-Mead

        # function evals   = 1618

        # data points      = 100

        # variables        = 7

        chi-square         = 113010.931

        reduced chi-square = 1215.17130

        Akaike info crit   = 717.006964

        Bayesian info crit = 735.243155

    ##  Warning: uncertainties could not be estimated:

    [[Variables]]

        C:        4891.30974 (init = 5000)

        I_inc:    1.96940567 (init = 0.4)

        sigma:    1.01340684 (init = 1)

        d_sigma:  0.08090573 (init = 0.07)

        phi:      0.09836989 (init = 0.2)

        kappa:    0.12588336 (init = 0.18)

        A:        9.30818331 (init = 10)

%% Cell type:code id: tags:

``` python

```

## load fitting data

def load_IQ_file(file_name):

    data = np.genfromtxt(path+filename,delimiter=',')

    data = np.genfromtxt(file_name,delimiter=',')

    QD_exp = data[:,0]

    IQ_exp = data[:,1]

    return QD_exp, IQ_exp

    minner = Minimizer(lmbda, params, fcn_args=(QD_exp, IQ_exp, np.arange(len(QD_exp))))

    result = minner.minimize('nelder')

    # write error report

    report_fit(result)

    # report_fit(result)

    print('Printing results to output.txt')

    with open('/src/output.txt', 'w') as f:

        for param in result.params.values():

            string = "%s:  %f +/- %f (init = %f)" % (param.name, param.value, param.stderr, param.init_value)

            print(string)

            f.write(string)

if __name__ == '__main__':

    print(sys.argv)

    print(len(sys.argv))

    if len(sys.argv)== 5:

    if len(sys.argv)== 10:

        filename = sys.argv[2]

        initial_guess = sys.argv[3:10]

        file_name = '/src/example_IQ.csv' # testing

        fit(initial_guess, file_name)

        # file_name = '/src/example_IQ.csv' # testing

        fit(initial_guess, filename)

    else:

        print("Usage: python this_code.py [-i] filename C I_inc sigma d_sigma phi kappa A")

        print("Usage: python3 /src/curvefit.py [-i] filename C I_inc sigma d_sigma phi kappa A")

        # python3 src/curvefit.py -i /src/example_IQ.csv 4100,4000,5000 0.5,0.1,2 0.975,0.9,1.1 0.07,0.05,0.1 0.08,0.05,0.5 0.15,0.1,0.18 10,0.5,20 

        sys.exit(1)

 No newline at end of file